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SECOND YEAR PROGRESS REPORT 


This report covers technical progress during the second year of the NASA Space Physics 
Theory contract “The Structure and Dynamics of the Solar Corona,” NAS5-96081, between 
NASA and Science Applications International Corporation, and covers the period July 16, 1997 to 
June 15, 1998. Under this contract SAIC, the University of California, Irvine (UCI), and the Jet 
Propulsion Laboratory (JPL), have conducted research into theoretical modeling of active regions, 
the solar corona, and the inner heliosphere, using the MHD model. During the period covered by 
this report we have published 17 articles in the scientific literature. These publications are listed in 
Section 4 of this report. In the Appendix we have attached reprints of selected articles. 

1. INTRODUCTION 

Our program centers around the theoretical modeling of active regions, the solar corona, and 
the inner heliosphere, using the MHD model. During the second year, we have placed increasing 
emphasis in advancing from qualitative comparisons of our results with observations to 
quantitative comparisons. Quantitative comparisons are important because they allow us to deduce 
solar parameters that are beyond our theoretical understanding at present (such as coronal heating). 

Our long-term goal is to develop a model that can be compared directly to observations with 
increasing fidelity. 

We have steadily increased the number of our collaborations with other solar physics groups. 
These have been largely driven by the success of our global coronal model. To date, we have 
compared our models with WIND and Ulysses spacecraft data, IPS measurements, radio 
emission, energetic particle events, coronal hole estimates, eclipse and coronagraph images, and a 
magnetic reconnection laboratory experiment. 

During the second year we have added the capability to follow the real-time evolution of the 
solar corona as it responds to changes in the photospheric magnetic field. We have used this 
capability to study the evolution of the corona during the rising phase of the new solar cycle (see 
Section 2.1.1). We have also made some important theoretical advances in the basic theory of the 
initiation of coronal mass ejections (CMEs). We have shown that emerging magnetic flux can be a 

very effective mechanism for the destabilization of coronal streamers and active regions (see 
Sections 2.2.1 and 2.3.2). 


2. ACHIEVEMENTS 

In this section we summarize the accomplishments made by our group during the second year 
of our Space Physics Theory Program contract. The descriptions are primarily intended to 
illustrate our principal results. A full account can be found in the referenced publications. 

2.1. Modeling the Large-Scale Structure of the Corona and Inner Heliosphere 

We have previously described an improved thermodynamic MHD model (in the First Year 
Progress Report) that has considerably extended the realism with which we can model the large- 
scale solar corona. In this section we describe the application of this improved model, as well as 
our polytropic model, to studies of the structure and dynamics of the solar corona. 
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2.7. j. Coronal Evolution in R esponse to Changes in Photnspheric Magnetic Fields 

Frequently we use our model to compute steady-state coronal solutions for a given 
distribution of photospheric magnetic field (e.g., as supplied by a synoptic map of the radial 
magnetic field from observations; see Secs. 2.1.2 and 2.1.3). This approach is limited to the study 
of the long-time properties of the solar corona. In reality, even if we neglect large-scale eruptions 
like coronal mass ejections, the corona is changing continuously, even during times of solar 
minimum, as evidenced by the recent Whole Sun Month campaign and as seen in SOHO 
observations. This changing structure is driven by changes in the photospheric magnetic field; 
active regions emerge and disperse continuously during the solar cycle. We have extended our 
model to incorporate the evolution of the photospheric magnetic field (a boundary condition for our 
model), so that we can now follow the evolution of the corona. This has given us the capability to 
study the long-term evolution of the corona (see below), the detailed evolution during a time period 
of interest (e.g., during Whole Sun Month), as well as the ability to study theoretically the coronal 
consequences of the emergence of magnetic flux (see Sec. 2.2.1). 

When we seek steady-state solutions (of Eqs. (1)— (9) in the First Year Progress Report), we 
set the tangential component of the electric field at the boundary, E fo , to zero. This keeps B ro (the 
radial magnetic field at the solar boundary) fixed in time. To make the flux evolve to match 
observed changes, it is necessary to specify a non-zero E/ 0 . In general, E ?0 can be expressed as 
V ? x ‘F r + V,0, where *F and 0 are arbitrary functions (of 0 and (p) and V, indicates tangential 
derivatives (in the 6-<j) plane at r = R s ). The potential 0 changes E/ 0 without changing the flux 
B ro , and can be used to control the transverse magnetic field (i.e., the shear and the normal electric 
current), whereas the potential *F changes the flux. We first consider the case 0 = 0. In this 
case, cV t 1 'F=dB r Jdt, which can be solved for ¥' for the flux change specified by dB r Jdt. 
Therefore, *Fis evaluated as new solar magnetic field measurements become available, specifying 
the evolution of E, 0 , which is used as a boundary condition for the MHD equations. 

Thus, rather than computing a sequence of steady-state solutions for each set of magnetic 
field boundary values, our time-dependent MHD model now represents the actual state of the 
corona corresponding to the evolving magnetic field measured on the surface of the Sun. 

It is important to note that a fundamental aspect of the state of the solar magnetic field is not 
provided by line-of-sight magnetograms. Since line-of-sight magnetograms do not provide 
information about the transverse component of the magnetic field, configurations with different 
levels of magnetic field twist cannot be distinguished (Mikic & Linker 1997). A vector 
magnetogram is required to uniquely specify the magnetic field (e.g., Mikic & McClymont 1994). 
In general, magnetic twist may change as a result of magnetic flux emergence or shearing flows at 
the photosphere. In principle, full-disk vector magnetograms can provide information about the 
evolution of the transverse magnetic field. 

We have used a sequence of synoptic Kitt Peak observations to study the evolution of the 
corona during the period Feb. 1, 1997-Mar. 18, 1998 (15 Carrington rotations). This time 
interval covers the beginning of the new solar cycle, as the Sun emerges from solar minimum, with 
the appearance of high-latitude active regions. To model the evolution over a time interval of over 
a year is computationally prohibitive at present. In order to study the quasi-static evolution of the 
corona, we changed the photospheric magnetic field at a rate that was enhanced by approximately 
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ten times compared to real time. (This simulation took ~ 65 hours of CPU time on a single 
processor of the Cray-T90 supercomputer. A real-time solution would take about 10 times longer. 
It may be possible to perform a real-time simulation once our code is ported to a massively parallel 
computer; see Sec. 2.4) This approximation makes it impossible to study the detailed evolution of 
individual events, though it is still meaningful to study the quasi-static evolution of the large-scale 
structure of the corona. For this simulation we used a “minimum shear” specification, <P= 0, 
since measurements of the transverse field do not exist. 

Figure 1 shows the evolution of the streamer structure, the coronal hole boundaries, and the 
heliospheric current sheet during this time period. Note the increase in complexity of the coronal 
magnetic field as the Sun emerges from solar minimum. Movies of the evolution of the corona 
were presented at the 1998 Spring AGU/SPD meeting in Boston (Tarditi, Linker, & Mikic 1998). 

2.7.2. Comparison with Eclipse. Observations 

We have continued our tradition of using our 3D MHD model to predict the state of the 
corona during forthcoming total solar eclipses. (Our eclipse comparisons can be seen on our Web 
page, http://iris023.saic. com:8000/corona/modeling.html.) We have recently made a new 
prediction for the eclipse of February 26, 1998, observed in the Caribbean. In Figure 2 we 
compare our prediction to an observation by HAO. 

We are planning to predict the state of the corona during the forthcoming total solar eclipse in 
August 1999, which will be seen in Central and Eastern Europe, the Middle-East, and Western 
Asia. This prediction will be our most challenging yet, since this eclipse will occur close to solar 
maximum, when the structure of the corona will be considerably more complicated than in 
previous eclipses that we have simulated. 

2.7.5. M HD Modeling for Whole Sun Month 

MHD modeling is particularly useful for studying the solar corona and solar wind when a 
coordinated set of observations is available, since it can help to synthesize different measurements 
into a coherent picture. The Whole Sun Month campaign (WSM; Aug. 10-Sep. 8, 1996) brought 
together a wide range of space and ground-based observations during solar minimum. Our 3D 
MHD model played a central role in the interpretation of these observations. Comparisons of our 
model with WSM observations are described by Linker etal. (1998), Breen etal. (1998), Gibson 
et al. (1998), and Posner etal. (1998). 

Figure 3 shows a comparison between synoptic images from our model and white-light 
observations from the Mauna Loa MK3 coronameter and the LASCO coronagraph aboard SOHO. 
The simulated and actual synoptic images show many similar overall features, indicating that the 
structure of the streamer belt during WSM has been captured by the MHD model. Both sets of 
images show a dark region breaking up the streamer belt near 270° longitude; this phenomenon is 
associated with the “elephant’s trunk,” an equatorial extension of the northern polar coronal hole. 
The elephant’s trunk was perhaps the most conspicuous coronal feature observed during WSM, 
and was apparent in several different wavelengths, including SOHO EIT images, NSOKP He 
10830 maps, and Yokhoh soft X-rays. It was most visible around August 26-27, 1996. Figure 4 
shows tracings of the magnetic field from the MHD model as they would appear on August 27, 
1996, with coronal holes (black) and closed-field regions (gray) mapped on the surface of the Sun. 
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Evolution of the Solar Corona During Feb. 1997-Mar. 1998 
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Figure 1. The changing structure of the solar corona during the period Feb. 1997-Mar. 1998 (Carrington rotations 
1919-1934), as illustrated by coronal hole maps (longitude vs. latitude, with gray/black indicating closed/open field regions), 
field line traces with the radial magnetic field shown on the surface of the Sun, polarization brightness, and the shape of the 
heliospheric current sheet. This time period represents the rising phase of the new solar cycle. The photospheric magnetic 
field was set as a time-dependent boundary condition on the 3D MHD simulation using Kitt Peak synoptic maps. 
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February 26, 1998 
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Figure 2. Comparison of an MHD prediction of the solar corona with a total solar eclipse observation. 
The prediction of the state of the corona was done using a 3D MHD simulation with boundary conditions 
on the magnetic field supplied by Kitt Peak observations collected during the period January 18 - 
February 12, 1998. The eclipse was visible in the Carribean. The eclipse image is courtesy of the the 
High Altitude Observatory, NCAR, Boulder, Colorado, USA. NCAR is sponsored by NSF. 
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MHD Modeling of the Solar Corona During Whole Sun Month 
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Figure 3. Comparison of simulated and actual white light synoptic images. Synoptic images are white light measurements 
that have been assembled throughout a rotation into a chart of brightness vs. latitude and longitude. The simulated and actual 
images show similar features, including the coronal hole that breaks up the streamer belt near 270 degrees. 
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Figure 4. Comparison of the MHD model with coronal holes seen in disk measurements on August 27, 1996. The "elephant’ 
trunk'' coronal hole (extending from the northern pole past the equator) can be seen in both the simulation and the data. 
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Figure 6. The solar wind speed measured at Ulysses and 
WIND mapped back to the Sun. The fast wind (red points) 
typically maps to deeper within coronal holes than the slow 
wind (blue points). 


Figure 5. Magnetic field lines traced back to the Sun from Ulysses 
(green) and WIND (blue), together with the heliospheric current 
sheet and Ulysses and WIND trajectories mapped back to 15 R s . 





For comparison, the NSOKP coronal hole boundaries and EIT images are also shown. It is 
apparent that the MHD model has captured the elephant’s trunk coronal hole, although the 
observations show the coronal hole extending to lower latitudes than predicted by the model. 

We also investigated the solar origins of features observed by the Ulysses and WIND 
spacecraft during WSM. Figure 5 shows tracings of field lines from Ulysses (green) and WIND 
(blue), as well as the heliospheric current sheet (HCS) predicted by the MHD model, and the 
Ulysses and WIND trajectories in the Sun’s rotating frame (mapped to the simulation domain). 
The model predicts HCS crossings by WIND (but not by Ulysses) during the WSM time period. 
WIND HCS crossings similar to those predicted were in fact observed. Figure 6 shows the 
observed solar wind speed, mapped back to its source at the Sun, with red points indicating the 
fastest observed speeds and blue the slowest. It is apparent that the slow wind maps back close to 
coronal hole boundaries, while the fast wind typically comes from deeper within coronal holes. 

This result was also seen during the Ulysses fast latitude scan (reported in the First Year Progress 
Report). 

2.7.4. Modeling the Solar Wind 

We have continued to perform extensive modeling in 2D (axisymmetric) geometry in order to 
self-consistently model the solar wind from its origins in the low corona to its expansion into 
interplanetary space. Our results show promising matches with generic in situ observations of the 
fast and slow solar wind, as well as the observed properties in the low corona. Full details were 
presented at the 1998 Spring AGU/SPD meeting in Boston (Lionello, Linker, & Mikic 1998). 

Our improved thermodynamic model significantly advances our capabilities to study the 
large-scale corona and inner heliosphere. In the new model, the high and low temperature regions 
in the corona are determined self-consistently by the open/closed field topology and the coronal 
heating profile. The closed-field regions naturally become hotter due to the insulating effect of 
parallel thermal conduction. Our 2D computations have shown that we can obtain the correct 
density contrast between coronal holes and streamers. As we move to a 3D model, we will 
investigate how the plasma density (through the polarization brightness) compares quantitatively to 
Mauna Loa and LASCO data for specific observations, and we will compare the solar wind 
parameters predicted by the computation (density, temperature, and velocity) with WIND and 
Ulysses data. 

2.7.5. Simulated Dis k and. Coronal Emission Imag es 

Our theory program has focused on the importance of developing diagnostics that are as close 
as possible to actual measurements, so that direct comparison between models and observations 
can be achieved whenever possible. The improved thermodynamic description in our MHD model 
opens up the possibility of modeling disk emission, just as we have previously done for 
polarization brightness. For example, as part of our comparisons for the Whole Sun Month 
campaign (Sec. 2.1.3), we compared coronal holes observed in EIT with the open-field regions 
predicted by the MHD model. While the dark emission features identified as coronal holes are 
clearly associated with magnetically open regions on the Sun, the observations do not directly 
measure field topology, and at times images in different wavelengths can yield different 
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interpretations of the coronal hole boundary. To make a quantitative comparison between coronal 
models and coronal hole observations requires reproducing the emission measurement itself. 

The emission lines observed by EIT arise from the excitation of iron ions, a trace species in 
the coronal plasma. The iron population in the corona is especially sensitive to temperature and has 
been modeled in the CHIANTI package (Dere et al. 1998) that is publicly available. Using the 
plasma temperature predicted by our improved MHD model, we can model the EIT emission for 
different lines. An example is shown in Figure 7. To generate this image, we used the polytropic 
MHD model (since our improved model is not operational in 3D yet), so that the temperature 
contrast in the model is not as large as expected. (This is a known deficiency of the polytropic 
model. We are merely using the polytropic model here to illustrate the technique.) The contrast 
observed in the synthesized emission (from the model) is produced by the contrast in plasma 
density (since the emission is proportional to the line-of-sight integral of the square of the electron 
density). We expect the agreement between our synthesized emission images and EIT 
observations to improve significantly when we incorporate the improved thermodynamics into our 
model, since this will allow the important effect of temperature on the iron population to be 
included. 

2.2. Coronal Mass Ejections 

2.2.1. Emerging Flux as a Trigger for CMFs 

Recently we have explored the effect of emerging flux on the stability and evolution of helmet 
streamers. In an influential paper, Feynman and Martin (1995) observed that CMEs (as deduced 
from disappearing filaments on the disk) are closely correlated with newly emerging magnetic flux 
in the neighborhood. Our emerging flux capability (described in Section 2.1.1) has allowed us to 
study the coronal effects of the emergence of subphotospheric magnetic flux. We have found that 
a sheared helmet streamer (i.e., with twisted magnetic field lines) erupts when magnetic flux with 
the opposite magnetic polarity emerges near the neutral line. This phenomenon has a threshold: 
when the amount of emerged flux is below a critical level, the helmet streamer is stable, and 
evolves quasi-statically; when a critical level of emerged flux is exceeded, the helmet streamer 
erupts violently, liberating almost all of its free magnetic energy, and ejecting a plasmoid outward 
with considerable kinetic energy (i.e., reminiscent of a fast CME). Note that it is necessary for the 
helmet streamer to be twisted for the eruption to occur — streamers with no twist do not erupt. 

Equilibria below the threshold naturally form a “filament” (or “flux rope”) above the neutral 
line. We have performed extensive simulations of this effect in 2D axisymmetric geometry (Mikic 
& Linker 1998), and we have begun to study the effect in 3D geometry. The (stable) “flux rope” 
(when the emerged flux is below the threshold) in 3D geometry is shown in Figure 8. The 
eruption in 3D geometry is shown in Figure 9. This is a very promising theoretical model for the 
initiation of CMEs by emerging flux. 

2.3. Three-Dimensional Active Region Modeling 

2.3.7. Interaction of Twisted Coronal Loops 

We have also modeled a more complicated situation in which two coronal loops interact 
during their emergence. We begin with a fully emerged twisted loop in an equilibrium state. A 
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Simulated EIT Emission 



Figure 7. (a) Disk emission at 195 A (Fe XII) from SOHO EIT on August 26, 1996. 
(b) Simulated disk emission using the polytropic MHD model (see text). 



3D MHD Simulation 


Figure 8. Formation of a stable flux rope 
within a sheared helmet streamer by the 
emergence of opposite polarity flux near the 
neutral line. Helical field lines like those 
shown here are believed to support 
prominences, which are often observed 
beneath streamers. 


Figure 9. Continued emergence of opposite polarity flux leads to the 
eruption of the flux rope and helmet streamer (i.e., a coronal mass 
ejection). No flux was emerged on the opposite side of the Sun, so this 
part of the streamer remains stable. 
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second loop, intended to model a rising flux tube emerging from below the photosphere, with 
different characteristics and physical orientation, emerges in the neighborhood. We have 
investigated several cases for the orientation of the second loop. In Figure 10 we show the 
evolution when the second loop emerges at an angle of 15° with respect to, and directly under, the 
first loop (with both loops having helicity of the same sign). We are continuing to investigate the 
details of the evolution and interaction of such loops. 

2-3.2. — Eruption of a n Active Region as a Result of Flux Emergence 

In Section 2.2.1 we showed that a twisted large-scale coronal arcade can erupt when flux of 
the opposite sign emerges near the neutral line. We have investigated the evolution of an active 
region arcade with respect to the same phenomenon, and we find that the arcade can erupt under 
similar circumstances. We first create a 3D twisted force-free field with a bipolar magnetic field 
distribution (by applying photospheric shear to a potential field), as shown in Fig. 1 la. We then 
emerge magnetic flux in such a way as to reduce the net flux in the active region (i.e., flux with 
opposite sign to the existing flux). When the amount of emerged flux is below a threshold, a 
stable “filament” is produced in the magnetic field above the neutral line (Fig. 1 lb). The twisted 
field lines in this 3D “filament” may be able to support a dense prominence. (Of course, this is not 
a filament in the traditional sense, since the model described here presently uses the zero-beta 
equations for simplicity. In the future, a more complete MHD model of this configuration may 
prove to be fruitful.) Note that the field line topology has not been prescribed: the formation of the 
“filament” is a natural consequence of the self-consistent evolution of the field, and has not been 
prescribed in any way. When the amount of emerged flux exceeds a critical threshold, the 
configuration erupts upward. We are continuing to study this configuration, which is a very 
promising candidate to explain the triggering of flares and possibly the eruption of filaments. 

2.4. Massively Parallel Computing 

Our spherical 3D MHD code was originally developed to take full advantage of vector 
computers (like the Cray-C90). In recent years, it has become clear that in order to perform more 
ambitious calculations (e.g., with increased spatial resolution and to simulate longer physical time 
periods), it has become necessary to use massively parallel computers (like the Cray-T3E). 
Unfortunately, it is not a trivial task to port an existing code to a massively parallel computer, since 

present compilers are not sophisticated enough to efficiently distribute the workload on multiple 
CPUs. 

To efficiently utilize a massively parallel computer it is necessary to parallelize the algorithm 
at a relatively low level. The de facto standards for parallel computing are High Performance 
Fortran (HPF) and Message Passing Interface (MPI). In order to decide which system was most 
efficient for our code, we parallelized a 3D potential field solver with HPF and MPI. This code is 
representative of our 3D MHD code, since both codes utilize a conjugate gradient solver for the 
inversion of matrices (where most of the CPU time is spent). 

Even though it took less time to implement HPF than MPI, our tests on the Cray-T3E show 
that the performance of MPI is superior to that of HPF. Ideally, the execution time (i.e., the 
average CPU time spent by each processor) should be inversely proportional to the number of 
processors. Overhead associated with communication between processors can cause deviation 
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Interaction of Two 
Coronal Loops 



Figure 10. The red loop is a twisted pre-existing coronal 
loop. The black loop emerges (with twist) underneath it 
and interacts with it. 


Initial Sheared Arcade Field Field with a "Flux Rope" after 

Flux Emergence 



Figure 11. (a) The initial sheared field in a 3D model of an active region. The magnetic field is in a force-free equilibrium, 
(b) After flux emerges near the neutral line (with opposite polarity to the existing flux), a stable equilibrium with a flux rope 
(black field lines) is formed. The arcade erupts if additional flux is emerged. 
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from this ideal “scaling.” Our results show that the scaling for MPI is closer to the ideal than for 
HPF. In fact, the CPU time for HPF “saturates” quickly as the number of CPUs is increased, and 
is not even able to match the execution time of the Cray-C90. (The CPU time for MPI also begins 
to saturate when the number of grid points per processor falls below ~ 300.) 

This evaluation has helped us to decide to use MPI to implement the parallel version of our 
3D MHD code. Even though implementing MPI is more complicated than implementing HPF, 
obtaining optimum efficiency is our most important concern. MPI appears to be the emerging 
standard for massively parallel computation, and is available on a wide variety of machines. 

2.5. Simulation of the MRX Laboratory Reconnection Experiment 

We have begun to apply innovative numerical techniques based on unstructured triangular 
grids to solar and astrophysical phenomena. We have successfully simulated a laboratory 
experiment that has been designed to study magnetic reconnection, especially as it relates to 
astrophysics. This technique is based on adaptive grids in which the mesh resolution is 
dynamically refined and coarsened to best represent the evolving solution. 

The MRX experiment at the Princeton Plasma Physics Laboratory (Yamada et al. 1997) has 
been constructed to systematically investigate the fundamental physics of magnetic reconnection in 
a well-controlled laboratory setting. The boundary conditions can be controlled externally, and the 
Lundquist number can be varied by controlling the plasma density and the magnetic field strength. 
Magnetic reconnection is induced by the forced merging of two toroidal plasmas. The plasmas are 
hot enough to be in an interesting parameter regime CIO 2 < S <: 10 3 ), but cold enough to allow 
the insertion of diagnostic probes to make detailed internal measurements of the magnetic field and 
plasma properties. This experiment, which is funded in part by NASA, can shed light on the 
fundamental processes involved in magnetic reconnection. 

We have used the nonlinear MHD code TRIM (Schnack, Lottati, Mikic, & Satyanarayana 
1998) to carry out direct numerical simulation of the MRX experiment. TRIM solves the time- 
dependent resistive MHD equations in axially symmetric geometry (although the plasma dynamics 
may be fully three-dimensional). TRIM uses an innovative algorithm based on an unstructured 
mesh of triangles. This allows the geometry of MRX to be accurately captured. TRIM also uses 
adaptive mesh refinement to dynamically refine and coarsen the mesh in the vicinity of evolving 
small scale structures, such as current sheets. 

The simulations show that the reconnection rate yas a function of the Lundquist number S 
scales as 5~«, with a ~ 0.3, over the range 10? < 5 < 10* implying that the reconnection rate 
is faster than the Sweet-Parker rate. We have used dynamically adaptive gridding to enhance the 
accuracy of these calculations. Our simulations show that a secondary magnetic island forms 
during the reconnection. A similar structure is seen in MRX under certain operating conditions. 
The secondary island is not observed in our calculations unless we use the enhanced resolution that 
is enabled by the adaptive grid. 
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immersed in a uniform-pressure plasma still remains an eqmlrbnum.) Srnce the magnenc 
field strength in the Gold-Hoyle equilibnum falls with the radial distance from the axis the 
plasma beta tncreases with the rad, us. reaching P = ! a. r - 10». We have chosen a finite 
beta to make the nonlinear evolution correspond more closely to the solar corona, w er 

the plasma B is small, but finite. • fi „u 

Since we expect the twist in the solar corona to be introduced by slow magnetic field 

footpoint motions in the photosphere, the appropriate initial state ought to be one m which 

the twist is only slightly larger than the critical value for linear instability. 

more footpoint shearing at z = ±L/2 and the applied surface flow considered ^ V - 0^ 

Eq (6)), because the time-scale evolution of the configuration (due to the kink) is smaller 
than Ihe one for the photospheric flow. We selected a iwis. of * = 3* m order do pro^ce 
a distinguishable linear phase of the instability, while at the same time keeping the excess 
twist (i.e., that above the stability threshold) small. In other simulations we have found that 
the nature of the nonlinear state does not seem to be sensitively dependent on the excess 
twist as long as it is not too large [26]. The loop length is set by the condition L 
giving a loop with aspect ratio L/a= 9.42 in this case. The (radial) AlWenhme is defin 

by r A = a/u^, where the Alfven speed on the axis is given by v A - , Um 

viscosity is used, corresponding to a viscous dissipation time r v = a /v - lUUr A . 

We study the ideal MHD evolution of this equilibrium (with r, = 0) in order to investigate 
whether the nonlinear evolution of the kink instability leads to the formation of curren 
sheets When strong gradients develop in the magnetic field during ideal MHD numer- 
ical simulations, it may be necessaiy to introduce plasma resistivity. In the case of the 
Gold-Hoyle field, as discussed below and as noted by Baty and Heyvaerts [25] the non- 
linear evolution of the kink instability does not introduce current sheets, so that it is not 
necessary to introduce resistivity into the calculation. This is in contrast to other equihbna 
that we have studied, of which the zero net-current equilibrium is a particular example 
for which we have found that the nonlinear evolution of the kink leads to the forma 
of current sheets [26], requiring the introduction of finite resistivity during the later s age 
of the calculation. We were thus able to perform the present calculation with the idea 
MHD model. (We note that a small amount of numerical resistivity is introduced during 
calculation by the upwind treatment of the advection, as described in Appendix BO 

We start the calculation at 1=0 with the m=0 equilibrium field given by 
Eqs (58)-(60), to which we add a small m = 1 perturbation with an amplitude v ~ 3 
1 J-4 w o ( T h e perturbation was chosen to be the eigenfunction corresponding to the most 
unstable linear mode in a periodic cylinder, modified suitably to have zero displacement at 
the axial boundaries, as required by line tying. Any small initial perturbation could 
been used without affecting the nonlinear results.) The equations were integrated for 500 r A 
requiring about three CPU hours on the Cray YMP/C-90 at NERSC. This code has also 

been implemented on the Cray T3D at CINECA in Bologna. 

The initial time step was chosen to be 0.1r A . The time step remained 0.1r A during the 
linear part of the run, decreasing to 0.05 t a during the initial phase of the nonlinear ^evolution 
as a result of the advective flow limit on the time step (Eq. 19) and increasing back to 0. 1 r A 
after saturation of the kink toward a new equilibrium. The advantage of usmg he semi- 
implicit scheme is illustrated by the fact that the wave Courant number (i.e., the rabo of the 
time step to the time step required by an explicit calculation) remains significantly larger 
than 1 during this calculation. Initially, the wave Courant number is 16, and it increase 
30 by the end of the calculation. 
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Thus we write 


v r;\,i = 


v 0 Xj 


Vz-Xj = 


1° 

if m ^ 1 

\~ iv e-,2j 

if m = 1 

f -V 6 -. 2 J 

if m = 0, 2 
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otherwise 
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(43) 


The first method to implement the boundary conditions, as explained in Subsection 3.1, 
can be used. In fact, the coupling of v r and v e in our scheme does not spoil the self-adjoint 
nature of the operator. 


3.3. The Induction Equation in Cylindrical Coordinates 
Let us consider how to solve the induction equation 

1 3 A 


T] 3 1 


= -V x V x A + S, 


(44) 


where the diffusive “curl-curl” operator is self-adjoint. We are not interested now in the 
ideal part of Eq. (1), which is treated separately with the predictor-corrector and merely 
adds a source term to the right-hand side of 
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\T] 
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(1 -co)AtV x Vx ) A (/2) . (45) 


The equation above is the finite-difference and self-adjoint representation of Eq. (44) when 
S = 0. It can be shown that to have stability for all At it must be \ < co < 1. In the code 
we choose co = which is also second-order accurate in time. The analytical form of the 
curl-curl operator, in cylindrical coordinates and after a Fourier transform in 6, is 
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Vz-ij - 1 = S z -ij. (36) 


Let us consider the boxed line in Eq. (35); it contains r*_ i at denominator which diverges 
for i = 2. The boxed derivatives in Eq. (32) produce that term. Using the results in Ap- 
pendix A we can rewrite the representation for m = 1 of the second diverging derivative in 
Eq. (32) as 


3 v r 
dr r 


— — + a\ + 0(r 2 ) = 

Yh- 2 


U;2 /y 2 ~ Y 
dY h ; 2 


V r; 2 = ao + <*irh;2 + 0(r 4 ). 


Now Y can be easily found, taking into account that 

2r h ;i = y 2 = dr h . 2 . 


(37) 

(38) 


(39) 


The first derivative can be treated in the same way. Hence, the boxed line of Eq. (35) when 
i — 2 and m = 1 becomes 


-\-ico] At (oi 2 j 6t] j) 



(40) 


When m = 2 we can write 


V r 


r 


dv r 

dr 


.dve _ .V0;2j — V0 ; \j 

dr { dr] 


The boxed line of Eq. (35) then becomes 


(41) 


a ./ , - \ V 0;2 ,j ~ V0.]j- 

-cO]At(a 2J +oi]j)r h -2 • (42) 

dr ] 

For m > 3 the above-cited line equals zero for i— 2. It is quite easy to verify that the 
equations for the components of v have been written in a self-adjoint form. Once we have 
fixed the boundary conditions we can deliver the equations to the CG solver. 

The field v may assume (in principle) arbitrary values V at z = =b L/2 and r = R, while 
for the singular boundary at r = 0 the values are dictated by conditions in Appendix A. 
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And, finally, for the z component we obtain 


l d dv z 

V • a'V\\ z — ra — - 

r dr dr 
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d dv z 

— a 

dz dz 


( 33 ) 


We choose as a natural grid for a the same as that used for po- This grid is shown in 
Fig. 4. When interpolated values are needed on other grids we put a A or a ' over the 
symbol to indicate respectively interpolation along r or z. We can now write the numerical 
representation of Eq. (30). We shall drop the temporal index for brevity. The right-hand 
side is indicated with just an Si, but its formula is analogous to the left-hand side, the only 
differences being the use of the old values of v and the factor — co 2 , instead of co\. The 
equation for the v r component is 
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The one for the vq component is 
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Note that when w > 0 the equations for u, and ve are coupled. For the v z component we 
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3. Get (M-SF^)/, i =2, N — 1. The only nonzero components are (M-\F (/?) )2 = 
A1 +; i Vo and (M • \F (/?) )/v-i = M + , N V L . 

4. Find boundary contributions to the right-hand side of (27), Sj — Si — (M • \F (/?) )/, i = 
2, N - 1. 

During the main iteration loop we calculate the “homogeneous part” of the solution, setting 
V = 0. At each iteration of the CG solver the matrix-vector product (M • \F)/ is evaluated 
in two steps: 

1. Set the boundary points on \F. Since V = 0, \F i = — C 0 ^ 2 , and N — —Cl'&n - l 

2. Get (M • \F);, i = 2,N-1. 

After the CG algorithm iterations, when the solution is retrieved, its values for i = 1 
and i = N are easily deduced from Eq. (28). 


3.2. Viscosity Equation in Cylindrical Coordinates 

In order to solve Eq. (4) the viscosity and semi-implicit operators must be inverted. If we 
write the finite difference equivalents of Eqs. ( 1 6)-( 1 7) as we have the numerical counterpart 
for the one-dimensional model in Eq. (24), we obtain matrices that are no longer self-adjoint 
as the analytical operators and require a slow algorithm for their inversion. However, it 
is possible to implement a general numerical scheme for inverting both operators. We 
shall represent them as symmetric and positive-definite matrices, using Subsection 3.1 as a 
guideline. We write 

(, pdV - Atoi dW • aV)v (,!+1) = [pdV + Atco 2 dV V • aV]v (,,) . (30) 


For the viscosity operator a is vp (we write p instead of po to avoid subscript overloading). 
Since we use a fully implicit advancement to ensure an efficient damping of the small 
unresolved scales, we take co\ = land C 02 = 0. For the semi-implicit algorithm a = C 2 Atpo 
and co\ — 1 , 02 = — 1 . Note that it is necessary to multiply both sides of the equation by the 
element of volume dV in order to obtain a self-adjoint matrix. This enables us to use the 
CG algorithm. The differential operator d W • QfV must be represented in a discrete form 
without spoiling the symmetry of the matrix. Thus, we first write this operator in cylindrical 
coordinates in the following form: for the r component it is 


V • afVv| r = 
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dr r dr dr r 
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im 
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For the 6 component it is 
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Note that in this form the matrix M is explicitly symmetric. Here we have dropped the 
temporal index and have indicated the known term just with Si. In computing (M • ^) z 
(where is either the guess solution or a temporary vector used by the CG algorithm), the 
index i may vary only in the range 2 < i < N — 1 . When we know how to fix the proper 
boundary conditions, without spoiling the self-adjoint nature of the operator, we shall be 
able to apply the CG algorithm. 

Boundary conditions are generally Dirichlet conditions (values specified at boundary 
points) or Neumann conditions (normal gradients at the boundaries). But even a combination 
of both is possible if we write them as 


+ Co ^2 = Vo forv = 0, 
&N i = V L for x = L. 


(28) 


Where Vo,l and C 0 ,l are constants whose values determine whether Dirichlet (C = 0) or 
Neumann (C = —1) conditions apply. A proof that the matrix M is positive definite, with 
constrains as in Eq. (28), is given in Appendix C. 

The way we implement the boundary conditions in our algorithm depends on how we in- 
tend to calculate the matrix-vector product at each iteration. We shall examine two methods: 
the first one modifies the matrix itself, the second one modifies the vector \E Z . 

To implement the first scheme we store the main diagonal A4o ; , and one offset diagonal 
A 4 + -j of the matrix. The diagonal values and the right-hand side S z must be modified 
according to Eq. (28) when i = 2 and when i — N — 1 as 


Ad();2 = A4 o;2 — C$AA + - 1 , 

S2 — S2 — Vo 2 W+;i; 

Mq\N-\ = A^0;/V-1 — ClA4+-n, 

§n~\ = S N - 1 — VlM+-n> (29) 


For 3 < i < N —2, §i = Sj , and Mo-j = Mq-j . Now only values of ^ with 2 < i < N — 1 
enter into the calculation of M • \I> and M is a symmetric positive definite matrix. After the 
CG algorithm iterations, when we have found the solution for the internal points, we set 
<I>i and $> N according to Eq. (28). 

However, the previous method might not work when we deal with 2D or 3D problems and 
more complicated operators such as “curl-curl.” Then it is possible that boundary conditions 
couple two different components of a vector field. In those cases, when we cannot write the 
modified diagonals for the self-adjoint matrix as in (29), we rely on the following method 
that has the advantage that we do not need to explicitly write the diagonals. 

We have implemented a subroutine to set the boundary points of the vector \I> according 
to Eq. (28). Another subroutine calculates (M • \I>)/, receiving an N component vector in 
input and yielding an N — 2 component vector in output. We split the calculation of the 
solution of Eq. (27) into two parts. First, we calculate the “inhomogeneous part,” fixing 
the boundary condition for the right-hand side. We take finite V terms and use as a 
work-array. The steps we perform are: 

1. Set^ (6) =0. 

2. Set the boundary points on '& {,>) . Since internal points are all zero we obtain 'b'/” = V 0 
and 1 ®'Sv ) = V L . 
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FIG. 7. Mesh used to represent Eq. (22). The physical region is between 0 and L. <I> lies on the mesh whose 
points are marked with x ; (3 lies on the one indicated with o’s. dx t and dx h .j are the distances between neighbor 
points for each mesh. 
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where the distance between two mesh points on the half-integer and on integer mesh are 
indicated with dxj and dxh-j, respectively. If we write Eq. (24) in matrix form as it stands, 
we find that the matrix A is not symmetric. Hence, we cannot apply the CG algorithm 
(actually the matrix is tridiagonal, and we might use a fast ad hoc direct solver for such 
cases. However it loses this property when we increase the number of dimensions). 

We know that for functions that are zero at the boundaries the following equality holds: 


[ XDp&dx = [ <$>DpXdx. 
Jo Jo 


(25) 


We can write the numerical representation of Eq. (25) as a product between matrices and 
vectors, 


X dxD^ $ = $ dxD^ X. (26) 

Here dx is a diagonal matrix whose elements are dx^j . With the discretization given by (24), 
we find that dxD^ is a self-adjoint matrix. When is positive, the matrix is also positive 
definite. 

Therefore, if we multiply both sides of Eq. (24) by dxh-j, we obtain 


(M • d>)/ = Xio ; /d>/ + .A4+ ; /d>/+i + — Si m , 

Mo-i — dxh-i + co At 
Pi 

M+-i = —co At - — ; 

dxi 


P> | ft- 1 

dxi dxi-\ 


( 27 ) 
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3. SELF-ADJOINT REPRESENTATION OF THE DIFFUSIVE 
AND SEMI-IMPLICIT TERMS 

The differential operators in the MHD equations have the property of being self-adjoint. 
Let us concentrate only on the diffusive terms that appear in Eq. (1) (resistive diffusion 
operator) and in Eq. (4) (semi-implicit and viscous). When these equations are advanced in 
time implicitly the problem requires solving the algebraic equation 

Ax = b, (21) 

where A is the coefficient matrix, x is the unknown vector, and b is the known term. In our 
case the dimension of the matrix A is 31 J x 31 J at worst, when the equations for the three 
vector components are coupled. 

We shall show that it is possible to write for all the above-cited operators a matrix A that 
is self-adjoint and positive-definite. This preserves an important property of the analytical 
equation and has also a desirable numerical advantage; we can apply the conjugate gradient 
(CG) algorithm to rapidly compute the solution x, instead of more complicated and general 
methods. The theory of the CG method is given in [16], and an application to a problem 
similar to ours is in [17]. Briefly, the CG method is an iterative algorithm to find the 
solution vector of the linear system (21) through successive approximations. It involves 
the matrix A only in the context of matrix-vector multiplication. Differently from other 
iterative methods, estimates of the largest and smallest eigenvalues of the iteration matrix 
are not needed. However rapid convergence occurs when the ratio between the maximum 
and minimum eigenvalues of A (known as the condition number) is small. Since our matrix 
is diagonally dominant, we apply diagonal preconditioning and obtain a matrix with a 
smaller condition number. The techniques described in this section make the code about 10 
times faster than its previous version in [15], which uses the biconjugate gradient method 

[17]. 


3.1. One-Dimensional Model 

We present now a discussion about how to implement a self-adjoint representation of a 
diffusive operator in one dimension. We consider the following diffusion equation: 


_ a / a$\ 

dt dx \ dx J 




( 22 ) 


x is assumed to vary between 0 and L . We want to solve the equation above numerically. 
First we fix two staggered meshes as in Fig. 7. <I> is defined on the half-integer mesh (marked 
with x ’s), while ft lies on the integer mesh (marked with o’s). The grid points need not be 
uniform. A general finite difference method for solving Eq. (22) is 

<b n + ] _ <b n 

=ft)D«-$" +1 +(1 (23) 

At 

Operators and vectors are in bold when we refer to them as a whole; when we consider 
their components we write them in normal type. The value of co can be any number between 
0 (fully explicit) and 1 (fully implicit), co = \ corresponds to a (centered) second-order 
accurate in At time discretization. We rewrite the previous equation in components as 
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y «+l _ y ** 

At 


Vypo +l/2 Vv n + 1 

P o +!/2 


(17) 


We have marked with P’s and C’s respectively the predictor and corrector steps in the 
equations for A, p, p , and v. Quantities marked with a * or a ** index are provisional 
values needed in the predictor-corrector schemes and for a fully implicit treatment of the 
viscous term in Eq. (17). 

The constant C 2 is the semi-implicit coefficient. The semi-implicit method is described in 
[13]. This method is unconditionally stable with respect to all magneto-acoustic and shear 
Alfven modes. Hence, accuracy becomes the most relevant consideration in the choice of 
the time step. Briefly, the method consists of adding to the original momentum equation (4) 
a linear term multiplied by a coefficient proportional to the time step, 


1 

Po 


V 


9 9 3v 

■ C At poV-—. 

ot 


(18) 


This removes the small time-step restriction originally introduced by the wave term. Invert- 
ing the linear operator above is much less complex and requires less computer memory than 
using a fully implicit scheme. The advective terms in (l)-(4) are formally only first-order 
accurate in At, while the wave-like terms are second-order accurate (centered). The use of 
the semi-implicit method for the wave terms leaves only the stability condition 


\(kV) m&x At\ < 1, 


(19) 


due to the explicit treatment of advection. The quantity k is the magnitude of the largest 
wave vector compatible with the grid size at the point (r, , z ; , 0 k ), 


k = yjkj + £ 2 + £ 2 



( 20 ) 


Note the presence of the factor “3” in the expression for k$, due to the dealiasing algorithm 
which restricts the largest poloidal mode to M/3. 

In order to address the stability limits imposed by the advective terms in Eqs. (l)-(4) 
and to give a heuristic justification of Eqs. (19), (20), we present a one-dimensional Von 
Neumann stability analysis of the advection part of the algorithm in Appendix B. 

A stability analysis of our algorithm indicates that the wave-like terms are stable for any 
choice of time step, and the advective terms are stable when Eq. (19) is satisfied. How- 
ever, we have recently found that the coupling of the leapfrog advance of the wave-like 
terms with a predictor-corrector for the advective terms may introduce numerical insta- 
bility. This instability does not develop when there is sufficient viscosity in the algorithm. 
The calculations we describe in Section 4 have sufficient viscosity to prevent this numer- 
ical instability from occurring. We have analyzed this coupling, and we have devised an 
algorithm that does not suffer from this instability [19]. The fully implicit differencing of 
diffusive terms in Eqs. (1) and (4) does not introduce any stability limitation in the time 
step. 
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divergence of a curl and the curl of a gradient vanish identically. The components of B and 
J are naturally specified on integer or half-integer meshes according to their definition. The 
components of the nonlinear terms in Eqs. (l)-(4) are evaluated on the same grid of the 
field components on the left-hand side, using simple averaging where necessary. 

Physical boundary conditions for A and v are specified at z = ±L /2 and r = R according 
to Eqs. (5)-(6). 

Boundary values for p and p are not required to be specified in our formulation, but they 
can be evaluated for diagnostic purposes using extrapolation. At r = 0 we apply geomet- 
ric boundary conditions as shown in Appendix A. It is not completely straightforward to 
implement such boundary conditions, as they tend to spoil the symmetry properties of the 
operators we have to invert. See Section 3 for further discussion. 


2.2. Temporal Approximation 

The right-hand sides of Eqs. ( 1 )— (4) have advective, dissipative, and wave-like terms that 
are treated using predictor-corrector, implicit, and semi-implicit methods. We introduce a 
leapfrog time discretization for the various fields, defining A (together with p and p) and 
V at staggered time intervals. The resulting algorithm is 


A* - A n ~ ]/2 
At 

pp*+ 1/2 _ p^n- 1/2 

At 


V « x B «-l/2 ; 

= v" x B* 

V x V x A" +1/2 


-rj- 


V x V x A"- 1/2 


■n- 


p* - p n ~ l/2 

At 

pn +\/2 _ pn ~ 1/2 

At 

p * _ p n~ 1/2 

At 

pn +\/2 _ n - 1/2 

At 


v — v" 

At 


At 


= -V • (p n ~ 1/2 \ n ), 

= -V • ( p*\ n ), 

= -V • (p n ~ l/2 \ n ), 

= -V • (pV) 

-(r - 1 )p n ~ l/2 v-v n , 
= - Y n . Vv", 

- — \ n • Vv* 


p 

C 


P 

C 

P 

C 

P 

c 


+ 


jn+ l/2 x g«+ 1/2 Vp n 


+ 1/2 


yi+1/2 


^n + l/2 


+ ■ 


V • C 2 A/ 2 Pq + 1/2 V (v** - v”) 


_£o 

A ^ +i/2 


(9) 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


( 16 ) 
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r = 0 r ( r — R 

FIG. 5. Mesh for v z , A : , B r , and J._. The square represents the physical domain in r and z. Mesh points are 
indicated with o. 

1 2 I - 1 



r = 0 


r = R 


FIG. 6. Mesh for B 0 ,v, and rj. The square represents the physical domain in r and z. Mesh points are indicated 
with •. 
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12 I -l 

o o o o J 



O O O O 1 

r = 0 t , r — R 

FIG. 3. Mesh for v r , A r , B : , and The square represents the physical domain in r and z. Mesh points are 
indicated with o. 


1 2 I 

X X X X X J 



FIG. 4. Mesh for v e , A 0 , J 0 , p, and p. The square represents the physical domain in r and z. Mesh points are 
indicated with x . 
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The reality of / requires that f m — f* m , where f* is the complex conjugate of f m . Moreover, 
/o and f m / 2 +i have zero imaginary parts. We apply (8) to the MHD equations (1)— (4), ob- 
taining a set of M nonlinear partial differential equations in the variables (r, z, f) describing 
the evolution of the Fourier components of A, v, p, and p. 

We evaluate the nonlinear terms in Eqs. ( 1 )— (4) with a fully dealiased pseudospectral 
algorithm, as described in [20] . The pseudospectral method consists of computing operations 
either in Fourier space or in real space, according to where it is more convenient. Thus 
multiplication is performed in real space to avoid convolution, and derivatives in 0 are 
evaluated in Fourier space. We use a fast Fourier transform to transform between the two 
representations. However, multiplication generates aliasing errors, due to quadratic and 
higher nonlinearities. Hence, we truncate the 0 -spectrum (dealiasing) and retain only two- 
thirds of available Fourier space. 

In order to simplify an implicit treatment, we assume that rj and v do not depend on 0 . 
This choice makes the implicit viscous and resistive operators linear in 0, and, consequently, 
poloidal modes decouple in Fourier space. 

We choose two staggered meshes for each nonperiodic direction, r and z. Beside being 
second-order accurate in calculating derivatives (when uniform meshes are specified), in this 
method boundary conditions are specified naturally: for the magnetic field only the normal 
component is specified, while the tangential one is computed, and for the electric field the 
tangential component is specified, while the normal component is computed. Moreover, the 
algorithm has the property that the longitudinal and transverse parts of vectors are effectively 
decoupled, so that initially vanishing longitudinal and transverse components will vanish 
all the time. A consequence of this is that V • B = 0. 

Current sheets may form during the nonlinear phase of instabilities in our simulations. 
We therefore allow the mesh points in the radial direction to have nonuniform spacing in 
order to have locally enhanced resolution in the proximity of the center of the loop. The 
axial mesh is normally (but not necessarily) uniformly spaced. Radial mesh points on the 
integer mesh are indicated with (r 7 , / = 1 , 1 — 1), where r\ = 0 and rj_\ = R. On the half- 
integer mesh we write (rn-j , / = 1,1). The relationship between the two set of mesh points 
is r ] ri — (r t + r 7 _ x ) / 2. We define also the finite increments (dr^j = r z — r z _i , i — 2 , 1 — 1) 
and (dr { = r/ z;z+ i — , / — 1 , 1 — 1). In the axial direction we define (zy , j — 1 , J — 1), with 

z\ = —L/2 and z j_\ = L/2, with analogous definitions for Zh-j, dzj , and dz h .j. Figures 3, 
4, 5, and 6 show how the dependent variables are defined on each mesh. 

Derivatives are defined on integer or half-integer meshes according to Table 1. With 
these definitions the gradient, divergence, and curl operator can be implemented so that the 


TABLE 1 

Differential Operators and Their Corresponding 
Finite Difference Representations 


Operator 

Integer mesh 

Half-integer mesh 

1 a 

rdO 

m 

‘T C,J 

m 

i — Di ; 

r h -a 

a 

E U - Ei-I ,j 

r+uj - Fi.j 

dr 

di'hj 

dr t 

a 

G,,j - G m _, 

H iJ+l ~ H UJ 

Yz 

dz h .j 

dz j 
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applying the condition at the boundaries z = db L /2, 


-rr = (V x B) f , (5) 

at 

v = V, (6) 


where the subscript t indicates the tangential component of vectors (the normal component 
of A is advanced as in Eq. (1)). Hence boundary conditions are specified only on the 
tangential electric field and normal magnetic field. Equations (5)-(6) are also valid at the 
radial boundary at r = R, where a conducting wall is present and V = 0 . The wall is placed 
far enough from the plasma not to affect the physics being studied. 

In order to translate from nondimensional to physical quantities we have to specify three 
normalization variables for length, magnetic field, and density. For example, when modeling 
coronal loops, we can set L 0 = 10 8 cm , B 0 — 10 G, and p Q = 10 -15 g cm -3 . Then fields and 
scalars in Eqs. ( 1 )— (4) can be measured in terms of 


V a = B 0 (4np 0 )- 1 ' 2 
J 0 = cB 0 (4nL 0 )~ l 
P 0 = B 2 0 (4 tt)- 1 

to =l 0 b~H^Pc ) 1/2 

r)o = (4jt) l/2 c~ 2 p~ l/2 B 0 L 0 
v 0 = B 0 L 0 (4np 0 )~ l/2 


G cm, 
cm s -1 , 
statamp cm -2 , 
dyne cm -2 , 
s, 
s, 

2 —1 

cm s . 


In terms of this normalization, we have: the Alfven velocity u A = V 0 ~ 9 x 10 9 cm s 1 , the 
Alfven time r A = L 0 / u A ~ 1 s, the mass scale M = p 0 L 3 = 1 x 10 9 g, etc. 


2.1. Spatial Approximation 

We use cylindrical coordinates (r, 0 , z), withO < r < R,0 < 0 <2ir , — L/2 < z < L/2, 
to model large aspect-ratio coronal loops. A sketch of the coordinate system is presented in 
Fig. 2. The 0 coordinate is periodic, so we introduce a discrete mesh 0j — 2n (j — 1 )/M, 
j = 1, 2, . . . , M, and write any field / as a finite Fourier series, 

M/2+1 

f(r,0j,z)- ^ f m (r, z)e lm8 >. ( 7 ) 

w = -M/2+1 

It is well known that the discrete Fourier series converges rapidly if the solution is smooth 
[18]. Furthermore, time advancement in Fourier space is facilitated because the poloidal 
(m) modes for linear operators decouple. Hence, implicit terms, which are present in the 
induction and momentum equations and must be inverted, will be represented as distinct 
small submatrices (one for each Fourier mode), instead of a single large matrix. 

The complex coefficients f m are given by 

1 M 

Mr, z) = -Y J f(r,e j ,z)e-‘ me >. 
j = i 


( 8 ) 
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photosphere 



FIG. 2. Coordinates and boundaries for modeling coronal loops. Field lines are anchored in the photosphere 
at z = ±L/2. The loops in the code are “straightened out,” since we ignore the curvature observed in real loops 
(see Fig. 1). 


The induction equation (1) allows us to select between ideal (rj = 0) or resistive MHD 
(: rj ^ 0). An ideal run is possible only in particular conditions, for example, to study the 
linear phase of an instability. In general, the grid resolution dictates the minimum values 
of r] and v that may be used. For example, for a 64 x 32 x 64 grid, we have found on a 
particular problem that they must be at least ^10 -3 , and sometimes ^10~ 2 , for the solution 
to be physically valid. 

When we discuss the zero-beta model, in which we assume that p — 0, we specify the 
density to be uniform and fixed in time, so that the mass continuity equation (2) is not 
solved. Similarly, we do not advance the energy equation (3) in the zero-beta model. Note 
that we neglect the influence of viscous and resistive heating, since we use an adiabatic 
energy equation. We plan to add the viscous and resistive heating terms, as well as thermal 
conduction, in future versions of the code. 

The viscosity in the momentum equation (4) is mainly used to damp short-wavelength 
modes in the calculation. In this term we have used p 0 = l/(27r) f pdO , instead of p, to 
allow the matrix inversion to proceed mode by mode. 

The equations describe the long-wavelength and long time-scale evolution of the corona, 
including magneto-acoustic waves, ideal and resistive instabilities, and resistive and viscous 
damping. However, particle acceleration and X-ray emission require kinetic models that 
are not part of the code. 

In a coronal loop, the dense photosphere anchors the footpoints of the magnetic field lines, 
so that they are dragged by applied surface flows V. This footpoint shearing is modeled by 
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a description of the implementation of a conjugate gradient algorithm for the inversion of 
the implicit spatial operators; in Section 4 we describe an application of the code, and in 
Section 5 we summarize our conclusions. 


2. COMPUTATIONAL MODEL 

Coronal loops consist of a hot, tenuous plasma embedded in a strong magnetic field. 
A most important feature of the loops is that their ends are firmly anchored in the dense 
photosphere. A sketch of a coronal loop is shown in Fig. 1 . 

The resistive MHD model is appropriate for our study of solar coronal plasmas. The 
MHD equations are written in cylindrical coordinates, neglecting for simplicity the curva- 
ture effect. Hence, loops in our analysis are “straightened out” as in Fig. 2. This is clearly 
an approximation and important effects are neglected in principle. This description is ex- 
pected to be appropriate when the aspect ratio (i.e., the ratio between the radial and the 
axial length scales) is large. We write the MHD equations in a convenient nondimensional 
form as 


dA 

~dt 


= vxB - rj'V x V x A, 


dp 

dt 


-V • (pv), 


— = -V • (pv) - (y - 1 )pV • v, 
ot 


3v 

Tt 


-v • Vv + 


JxB 

P 



V • vpoVu 
Po 


( 1 ) 

(2) 

( 3 ) 

(4) 


where A is the vector potential of the magnetic field B = VxA,J = VxBis the current 
density, v is the velocity, p the pressure, p the mass density, r] the resistivity, and v the 
viscosity. 



FIG. 1 . A schematic representation of the magnetic field of a loop in the solar corona. Note that all the field 
lines are anchored in the photosphere at both ends. In our code we neglect the curvature and the loop appears 
“straightened out” (see Fig. 2). 
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its structure and evolution. Loops of tenuous magnetized plasma are observed frequently in 
the corona. The ends of these coronal loops are anchored in the dense photosphere below, 
a situation that has been referred to as line tying. The slow motions in the photosphere 
drive the footpoints of the loops, which evolve through series of equilibria. Rapid evolution 
may occur when an unstable equilibrium is reached, possibly leading to the development 
of current sheets, magnetic reconnection, and the release of magnetic energy that may heat 
the corona. If such “disruptions” are sufficiently impulsive, they may be identified with 
solar flares. Recent observational and theoretical results on flares and coronal heating can 
be found in [1-4]. 

A complete description of these processes requires a three-dimensional model that in- 
cludes the slow, long-wavelength evolution prior to disruption, as well as the rapid short- 
wavelength evolution in the nonlinear phase. The resistive magnetohydrodynamic (MHD) 
model is appropriate to describe much of the physics associated with these phenomena (ex- 
cept in places where the gradient scale-length is smaller than the gyroradius and a kinetic 
treatment must be adopted). For simplicity we will restrict our attention to geometries that 
are best described in a cylindrical coordinate system. When modeling coronal loops, we 
will therefore neglect the important effect of loop curvature [5], studying instead straight 
flux tubes as an approximation to large-aspect-ratio coronal loops. 

Although the linear stability properties of cylindrical flux tubes have been studied ana- 
lytically [6-9], a description of the nonlinear evolution requires a computational approach. 
Several cylindrical MHD codes, with axially periodic boundary conditions, have been used 
to model laboratory plasmas [10-14]. However, we cannot use such codes for our studies 
because we need to impose line-tied boundary conditions to properly model coronal loops, 
and a new algorithm must be developed for this purpose. 

The goal of the present paper is to describe a fast, accurate, and reliable algorithm for the 
advancement of the full resistive and viscous MHD equations in cylindrical geometry which 
allows for the specification of driving photospheric motions at the magnetic footpoints. The 
code is an improved version of the algorithm employed in [15]. Quantities are evaluated on 
grids: the azimuthal variation (0) is represented using Fourier series, with pseudospectral 
calculation of derivatives; the r and z coordinates are discretized on staggered meshes, 
which allows us to define a curl operator whose divergence vanishes identically. A leapfrog 
scheme is used for the time advancement of the wave terms. We employ a semi-implicit 
operator in the momentum equation, following the method described in [13], while treating 
advection with a predictor-corrector scheme. The semi-implicit scheme allows us to set the 
time step through considerations of accuracy rather than stability of the algorithm and leads 
to a substantial saving of CPU time, compared to a fully explicit algorithm. 

The resistive and viscous diffusion terms are advanced implicitly. The resulting implicit 
equations and the semi-implicit operator are inverted using a preconditioned conjugate 
gradient method [16,17]. We have attempted to preserve many of the analytical properties 
of the MHD equations in the discretized equations. In particular, we have taken special care 
to preserve the self-adjointness of spatial difference operators. Since the proper differencing 
of a self-adjoint operator results in a symmetric matrix, we are therefore able to use the 
efficient methods that exist for inverting symmetric matrices. As an illustration of the 
properties of our algorithm, we describe its application to the linear and nonlinear evolution 
of a kink instability in a twisted flux tube. 

The paper is organized as follows: in Section 2 we describe the MHD equations and the 
spatial and temporal approximations employed to advance them in time; Section 3 contains 
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We describe a three-dimensional algorithm for the advancement of the resistive 
MHD equations in cylindrical geometry with line-tied boundary conditions. This 
code has been developed to simulate the behavior of solar coronal plasmas. A finite- 
difference discretization is used for the radial and axial coordinates; a pseudospectral 
method is used for the azimuthal coordinate. The dependent variables are defined 
on finite-difference meshes that are staggered with respect to each other to facili- 
tate the application of boundary conditions. The time-advance algorithm features a 
semi-implicit leapfrog scheme for the wave terms, a predictor-corrector treatment 
of advection, and an implicit advance of the resistive and viscous diffusion terms. 
The semi-implicit and implicit operators are inverted using a preconditioned con- 
jugate gradient method. Special care is taken in maintaining the self-adjointness of 
the discretized operators, so that a fast inversion algorithm applicable to symmet- 
ric matrices can be used. By way of illustration, we describe the application of the 
code to the linear and nonlinear evolution of a kink instability in a twisted flux 

tube. © 1998 Academic Press 

Key Words: partial differential equations; initial value and time-dependent initial- 
boundary value problems; numerical linear algebra; iterative methods for linear sys- 
tems; astronomy and astrophysics; hydrodynamic and hydromagnetic problems; fluid 
mechanics; magnetohydrodynamics and electrohydrodynamics. 


1. INTRODUCTION 

The solar corona abounds with interesting phenomena of controversial physical interpre- 
tation. Although it is not understood why the corona is so hot (around 10 6 K) and what causes 
flares to occur, it is believed that magnetic reconnection plays a crucial role in determining 
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In Figs. 8 and 9 we plot the magnetic energy and the kinetic energy in various modes as a 
function of time. Initially, the m = 0 mode shows the relaxation of the analytic equilibrium to 
the mesh (since it is not a perfect equilibrium of the discretized equations). The m — 1 mode 
grows exponentially in time, with a growth rate yx a = 0.022. The higher-ra modes show 
growth associated with the coupling to the m = 1 mode. Beginning at t ~ 200ta, when the 
m = 1 mode reaches a significant amplitude, there is a nonlinear interaction during which 
the higher-m modes become sizable. This phase corresponds to the observed kinking of 
the axis of the flux tube. Eventually, the kink mode appears to saturate, indicating that the 
kinked flux tube is settling toward a new equilibrium. 

The linear growth rate of the m = 1 mode at <I> = 3n is lower than previous estimates 
because of the effect of finite beta. For the case with zero beta, the growth rate has been 
estimated previously as yiA = 0.034 by Mikic et al. [15], yxp, = 0.027 by Foote and Craig 
[23], and yx& = 0.037 by Baty and Heyvaerts [25]. Apparently, even though the plasma 
beta is small on the axis, the growth rate is changed significantly by the plasma pres- 
sure. This is because the magnetic field strength falls far from the axis in this equilib- 
rium, so that even a small pressure can affect the kinking motion of the flux tube. Indeed, 
when we repeated the calculation with the zero-beta model (i.e., with p Q — 0 and a con- 
stant density), we found the linear growth rate of the m = 1 mode to be yx a = 0.038, in 
good agreement with previous zero-beta results. (The growth rate determined by Foote 
and Craig is only intended to be a rough approximation for this equilibrium near the 
marginal stability point [23].) The finite pressure leads to a reduction of the growth rate, 
apparently due to the fact that beta is greater than one at a large radius, as described 
above. The finite-beta case is a more realistic representation of the solar corona than the 
force-free case (with = 0), in which the flux tube kinking in the weak-field region is not 
impeded. 

Figures 8 and 9 show that the kinked flux tube appears to settle to a new equilib- 
rium state. This state does not appear to have any current sheets; the magnetic field re- 
mains smooth and free of discontinuities. In Fig. 10 we show the evolution of the total 



FIG. 8. Magnetic energy in various Fourier modes as a function of time for the nonlinear kink. The energy in 
normalized by the factor E 0 = B*a 3 / (8jt) . 
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FIG. 9. Kinetic energy in various Fourier modes as a function of time for the nonlinear kink. The energy in 
normalized by the factor E 0 = Bla 3 /(%n). 


magnetic, kinetic, and thermal energies (defined by W = f [B 2 /Stt] dV , K = f ^ pv 2 dV , 
and E = f [p/(y — 1 )]dV, respectively). Note that as the flux tube kinks, the magnetic 
energy is converted into kinetic energy and, finally, into thermal energy. The kinked flux 
tube approaches an equilibrium that has smaller magnetic energy than the initial state. 

The large-scale kinking of the flux tube is best illustrated by traces of the magnetic field 
lines. In Fig. 11 we show traces at four instants of time. At t = 100r A , during the linear 
stage, the kink is barely perceptible in the field line plot. At t = 250r A the kinking pattern 
is clearly visible. The traces at t — 400r A and t — 500t a show that the kink is saturating to 
a new equilibrium state. 
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FIG. 10. The total magnetic, kinetic, and thermal energies (indicated respectively with W, K , and E) as 
functions of time for the nonlinear kink. The energy in normalized by the factor E 0 = B 2 o a 3 /(Sn). 
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FIG. 11. Field line plots at t = 100t a , t = 250t a , t = 400t a , and t = 500r A . Field lines start from the bottom 
of the loop from a circle of radius a. Initially the kink pattern is barely visible, but when the instability saturates 
to a new equilibrium (third and fourth panel), the center of the loop has moved outward to about r = 4a. 


5. CONCLUSIONS 

We have presented a fast and accurate algorithm for the solution of the full resistive and 
viscous MHD equations in cylindrical coordinates in the presence of line-tied boundary 
conditions. The computer code based on this algorithm has been applied to the study of 
solar coronal flux tubes. In particular, the techniques are suited to the simulation of flux 
tubes whose footpoints are driven by slow photospheric motions. 

The algorithm is implemented using finite differences in two dimensions, with pseu- 
dospectral derivatives along the third (periodic) dimension. The use of staggered finite- 
difference meshes preserves the solenoidal nature of the magnetic field and leads to a natural 
specification of boundary conditions on the tangential electric field and the normal magnetic 
field. Time advancement of the wave-like terms is performed with a leapfrog scheme. A 
semi-implicit operator is used in the momentum equation to give unconditional stability 
to wave-like terms. Advective terms are advanced using a predictor-corrector scheme, and 
therefore limit the time step by a Courant condition based on the flow speed. This allows 
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us to use significantly larger time steps than those achievable by a fully explicit algorithm. 
The viscous and resistive terms are discretized using a fully implicit time advance. The 
semi-implicit, viscous, and resistive operators are inverted using a preconditioned conju- 
gate gradient method. Special care has been taken to maintain the self-adjointness of the 
discretized operators, so that a fast inversion algorithm applicable to symmetric matrices 
can be used. 

To illustrate the application of the code, we have presented the nonlinear evolution of 
the ideal kink instability in the Gold-Hoyle uniform-twist field. Our results show that it 
is possible to follow the linear and nonlinear evolution of the kink instability. In the case 
of the Gold-Hoyle equilibrium, it appears that the kink instability saturates nonlinearly as 
the flux tube evolves to a new kinked equilibrium without the formation of current sheets. 
This result is in agreement with the results of Baty and Heyvaerts [25]. In contrast, Craig 
and Sneyd [24] concluded that the kink instability in the Gold-Hoyle field causes current 
sheets to form, a conclusion that is based on a calculation on a Lagrangian mesh whose 
accuracy is impaired when the mesh becomes significantly distorted by the finite-amplitude 
kink displacement. The evolution observed here for the Gold-Hoyle equilibrium contrasts 
sharply with the nonlinear evolution of the kink mode in a tokamak in which the nonlinear 
evolution causes current sheets (i.e., true discontinuities in the magnetic field) to form [27], 
a difference that has been attributed to the effect of line tying in the case of the coalescence 
instability by Longcope and Strauss [28]. In our case, it was thus possible to study the ideal 
MHD evolution. In general, instabilities can introduce current sheets, in which case it is 
necessary to study the resistive evolution. Equilibria in which the kink instability creates 
current sheets are discussed in [26, 29]. The role of a resonant surface in the formation of 
current sheets as a result of the nonlinear evolution of kink instabilities has been addressed 
previously [9, 25]. 

Therefore, the kink instability in the Gold-Hoyle equilibrium is not likely to play an 
important role in the solar corona, since it does not appear to cause significant heating or to 
lead to impulsive motions. On the other hand, other equilibria, in particular those in which 
the nonlinear evolution causes current sheets to form, leading to significant plasma heating, 
magnetic reconnection, and particle acceleration, are likely to be of interest in understanding 
coronal phenomena. Numerical algorithms and codes such are the one detailed here will be 
an important tool in this endeavor. 

The code has also been used elsewhere [29, 26] to study the nonlinear evolution of 
instabilities in more realistic equilibria that are intended to model coronal loops formed by 
the twisting of uniform ambient fields and from the emergence of magnetic flux tubes from 
the photosphere. In these cases we have modeled the formation of current sheets, magnetic 
reconnection, and fast energy release. 


APPENDIX A: FOURIER COEFFICIENTS IN POLAR COORDINATES 

Let us consider a scalar function F(x, y). We assume it is regular near the origin and we 
expand it in Taylor series 
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F(x, y) = F 0 +x 
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Since x = r cos 0 and y = r sin 6, We can rewrite Eq. (62) as 
F(x, y) = F 0 +r 
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(63) 


Hence the r n term is a certain combination of exponential functions e im0 , with —n<m<n. 
Since p 7 " 70 may appear only when n > \m\ when we calculate the Fourier series of F, we 
obtain 

oo 

F (m \r) = (64) 

«=|w| 

Thus we have F (w) (r) = 0(r |m| ) for small r. 

Let us considernow only one term of the Fourier series a (w) (x, y) — F (m) (r)p 7m6> . Notwith- 
standing r is always defined to be greater than zero, we notice that in an algebraic point of 
view we are allowed to write 


r — r, 

6 — > 0 + 7T. 


In this case x and y do not change and so 

F (m) e im \r) = F {m \-r)e ime (-\) m . 


(65) 

( 66 ) 


(67) 


Let us expand both members of Eq. (67) in Taylor series around r — 0, obtaining 

CO 

Y, 4”V(l-(-l) n+m )=0. (68) 

This imposes the following condition on the Taylor series terms of the rath Fourier coeffi- 
cient of F(x, y). 
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(m) = 0 


(n+m = 2k + 1 
\ A: = 0, 1,2, .... 


(69) 


This means that the Taylor series of an even coefficient has only even terms and, vice versa, 
if m is odd only odd terms are found. 

Let us examine now a vector U = (U x , U y ), where the vector components U x (x, y) and 
U y (x , y) are scalar functions with the same properties of F(x, y). The components of U in 
polar coordinates are 


U r = U x cos 0 + U y sin 6, 
Uq — —U x sin 6 + U y cos 0. 


(70) 

(71) 
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From this follows that the Fourier series coefficients are 


Of = \ [uf 0 + uf +n - iOf l) + iuf +,) ) , (72) 

Of = \ [iuf~ X) - iuf +l) + uf~ X) + uf +]) ^ . ( 73 ) 

Thus Of and Of are 0(r mm when m is even only odd terms of the Taylor 

series are found and vice versa. Note that for r =0 and m > 1 the following equality holds: 

Of — — iOf. (74) 


In cylindrical coordinates the third component U z behaves as a scalar function. 


APPENDIX B: STABILITY OF PREDICTOR-CORRECTOR 
ADVECTION ALGORITHMS 


The typical advection equation in one dimension is 

9/9/ n 

b v — = 0. 

3 1 dx 


To solve the equation above we employ the predictor-corrector algorithm 
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where 0 < a < 1 . For centered differences 
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The scheme above is first-order accurate in time. In order to perform a Von Neumann 
stability analysis we suppose that a local solution behaves like / (j Ax , t n ) = z n exp (/ kj Ax) 
and we assume for simplicity that v > 0. The amplification factor z ( k ) must have modulus 
less than 1 for stability (see [30] for a more complete discussion of the method). Substituting 
/ into Eq. (76) we obtain 


vAt 


Ax 


vAt 


1-— Q [a — Q-l 


Ax 
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(79) 
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The case of centered differences and no predictor-corrector (a — 0) yields |z | always greater 
than one and is unconditionally unstable. If we introduce upwind differences then we have 
\z\ < 1 when 


At 


v < 1, 

Ax 


(80) 


the so-called Courant condition. With the predictor-corrector and centered differences we 
obtain 


2a — 1 / vAt \ 2 

a 2 ~ \ Ax ) 


( 81 ) 


for stability. For a fully advanced corrector (o' — 1) we find again the Courant condition. 
Using these methods introduces a numerical viscosity term 


V n 


9V 

dx 2 


(82) 


into Eq. (75), which is useful to damp small unresolved scales. Let us write z = exp(— ico r At 
+ y At) and then find y from Eq. (79), limiting ourselves to the case k Ax <C 1 . From Eq. (82) 
it follows that the numerical viscosity coefficient is v n = —y / k 2 . With upwind differences 
and a = 0 its value is 


vAx 

~ 2 ~ 



(83) 


and for a = 1 and centered differences 


The situation in the code is complicated, with respect to this simple example, by the presence 
of nonuniform three-dimensional meshes in a non-Cartesian frame of reference. Further- 
more, the conditions above are only necessary and not sufficient for stability. Fully advanced 
predictor-corrector is used to stabilize advection in the periodic direction 0, since we cannot 
upwind 0 -derivative. We normally combine this method with upwind differences in r and 
z, originating the stability condition showed in Eq. (19) since both must obey the Courant 
condition. 


APPENDIX C: POSITIVE DEFINITENESS OF A 
SYMMETRIC TRIDIAGONAL MATRIX 

Let us consider a symmetric tridiagonal matrix A of the form: 


c 2 + bi+b 2 + C 0 b i, 

Aij = ^ Ci + bt-i -\-bi, 

. Cv-i + b ^~ 2 + bfli-i + C^bu- 1 , 


Au + 1 =bi 2 < i < N — 2 
A U -i = bi-\ 3 < i < N — 1, 


i = 2, 

3 < i < N - 2, 
i = N — 1, 


(85) 
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where c, and bi are positive and Co and Cl may be either 1 or — 1 . A matrix A is positive 
definite if and only if 

x • A • x > 0 Vjc ^ 0. (86) 

The condition above becomes for our matrix 

N-\ N - 1 

Y,xf(« +bj +bi- 1) + YxjXi-ibj-i 

i = 2 /= 3 

N—2 

+ y ^ -U-U+i > 0 

i=2 

We can rewrite it as 
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x z 2 Ci + yy Z?/ (x/ — x /+ i) 2 + x 2 £i (1 + Co) + (1 + Cl) > 0, (89) 

/ =2 z =2 

that is manifestly true. 
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A new finite-volume algorithm for the solution of the time-dependent, nonideal 
magnetohydrodynamic (MHD) equations in cylindrical (r, 0, z) geometry is pre- 
sented. The boundary geometry is assumed to be axially symmetric, but it can have 
arbitrary shape and connectivity in the poloidal (r, z) plane. The dynamics of the 
fluid is fully three-dimensional. A two-dimensional, unstructured, adaptive grid of 
triangles is used to describe the poloidal geometry. A pseudospectral algorithm with 
fast Fourier transforms is used for the periodic toroidal (0) direction. The grid can be 
dynamically refined or coarsened by adding or deleting points to adapt to evolving 
fine-scale structures in the solution. The algorithm exactly conserves total mass, mo- 
mentum, energy, and magnetic flux, and identically preserves the solenoidal prop- 
erties of the magnetic field and the current density. Examples of the application 
of the algorithm to two-dimensional hydrodynamic and MHD shocks, the linear 
growth of a resistive tearing instability in a tokamak, and the linear growth and 
nonlinear saturation of three-dimensional kink instabilities in toroidal geometry are 

given. :;e) 1998 Academic Press 


1. INTRODUCTION 

Over the past two decades, large scale numerical simulation has played an important role 
in fusion plasma research. Applications of these techniques to fluid plasma models have 
led to an interpretation of sawtooth [1] and fishbone oscillations [2 — 4] in tokamaks, the 
tokamak major disruption [5], the tilting mode in field-reversed configurations [6], and to 
a fundamental understanding of the reversed-field pinch dynamo [7]. These calculations 
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were generally performed on spatial grids aligned with fixed coordinate directions. As a 
result, these calculations were often carried out in generic toroidal or cylindrical geometry. 

Many important engineering and theoretical issues are affected by the details of the 
geometry. For example, the poloidal plasma shape can greatly influence the linear stability 
properties of a fusion plasma, and these linear calculations are now routinely performed with 
the actual poloidal plasma geometry accurately represented. This is often accomplished by 
employing a coordinate system based on the magnetic field lines, whose geometry is fixed 
throughout the calculation. The resulting metric makes the fluid equations quite complicated, 
but it allows the coordinate system to naturally fit the plasma shape. 

While coordinate systems based on magnetic fields have proven useful in linear calcula- 
tions that principally determine eigenvalues and eigenvectors, they have several disadvan- 
tages in fully nonlinear simulations because of the dynamical nature of the magnetic field. 
Primary among these is the nonuniqueness of the magnetic topology when finite resistivity 
is included in the model. These coordinate systems also may become singular at magnetic 
separatrices, which are essential features of modem tokamaks. 

It is thus desirable for fully nonlinear simulations to employ a spatial representation 
that can readily conform to the geometric details of the plasma and its surroundings and 
is independent of the magnetic structure. For accuracy, this representation should also be 
capable of conforming to the dynamical evolution of short spatial scale structures, such 
as current filaments and density gradients that may appear spontaneously and require finer 
spatial resolution than the surrounding environment. 

One candidate for a spatial representation with these features is an unstructured , adaptive 
mesh. In such a mesh the mesh points are not constrained to lie along constant coordinate 
directions. Instead, mesh points are placed on the boundary to conform with the actual 
geometry of the problem and are distributed in space to maximize the accuracy of the 
calculation. Thus placed, the points are connected with line elements that form the edges of 
triangles (see Fig. 1). These triangles are the Eulerian control volumes that form the basis 
for the finite representation of the appropriate fluid equations. In the logical data structure 
that describes the mesh, mesh points (and associated triangles) can easily be added or 
deleted dynamically, based on predefined accuracy criteria. The spatial representation can 
thus adapt to evolving spatial structures without the mesh distortion problems associated 



FIG. 1. Triangle, edge, and vertex mesh elements. 



UNSTRUCTURED GRID MAGNETOHYDRODYNAMICS 


73 


with Lagrangian formulations. Also, since the placement of the mesh points is not tied to a 
magnetic field, separatrices and X-type neutral points can be readily resolved. 

Techniques based on unstructured, adaptive meshes have come to maturity in computa- 
tional fluid dynamics (CFD), where quantitative predictions in real geometry have become 
essential in the design of aircraft and gas turbine engines [8]. These methods are generally 
based the solution of a Riemann problem at each triangle interface (edge) to determine the 
fluxes of energy, mass, and momentum [9]. The simplest extension of the hydrodynamic 
model that is appropriate for the description of magnetic fusion plasmas is magnetohydro- 
dynamics (MHD). 

In this paper we describe an extension of these spatial gridding techniques to an MHD 
model suitable for the description of the fully three-dimensional dynamics of plasmas with 
axially symmetric geometry, such as toroidal fusion devices. Since the dominant MHD 
modes in this geometry have relatively long toroidal wavelength, the toroidal coordinate is 
approximated with finite Fourier series. The unstructured, triangular mesh is used to describe 
the details of the poloidal geometry. The hydrodynamic variables are treated in a manner 
analogous to that used in CFD. These quantities (mass, momentum, and energy) are volume 
based densities that satisfy scalar or vector conservation laws. The electromagnetic variables 
(the magnetic flux density B and the electric current density J) are area-based densities that 
satisfy pseudo-vector conservation laws and have no counterpart in fluid dynamics. These 
variables are constrained to remain solenoidal. These quantities are represented on the 
triangular mesh in a manner that is an extension of that used on rectangular, structured 
meshes. 

In this work we have chosen to solve the primitive (instead of reduced [10-12]) MHD 
equations in order to make the resulting codes and techniques more generally applicable to 
problems beyond the narrow scope of tokamak plasmas. The temporal stiffness problems 
inherent in this description of tokamak dynamics that motivate the reduced MHD model 
are addressed here with the semi-implicit method of time integration [13-16]. We remark 
that, while the present work deals strictly with the MHD equations, other volume-based 
fluid descriptions, such as diffusive transport, could easily be adapted to these techniques 
and coupled with the description of the electromagnetic field presented here. 

We emphasize that a primary goal of this work is to develop an algorithm for the descrip- 
tion of slow, nonlinear, long wavelength motions in toroidal geometry with arbitrary poloidal 
shape. We have therefore used several low-order approximations that may be inappropriate 
for problems in which the highly accurate description of strong shocks is required. The 
solution of such problems will require extensions of the work described here. Nonetheless, 
the present algorithm can reasonably describe shock formation and propagation at relatively 
low Mach number. 

This paper is organized as follows. In Section 2 we discuss the properties of structured 
and unstructured meshes, and the data structures useful for describing them. Issues related 
to the triangulation of an arbitrary set of points in a plane are also discussed. In Section 3 we 
derive a finite volume approximation to the resistive MHD equations suitable for use on an 
unstructured, triangular mesh in toroidal geometry. Boundary conditions are discussed here. 
The specific MHD model and its implementation on the unstructured mesh are discussed 
in Section 4. In Section 5 we discuss methods of time integration and describe our imple- 
mentation of semi-implicit and fully implicit algorithms. Examples of the application of the 
method are given in Section 6. Included are standard, two-dimensional hydrodynamic and 
MHD shock problems, as well as applications of the method to the stability and nonlinear 
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evolution of toroidal fusion plasmas in three dimensions. The summary and discussion are 
given in Section 7. 


2. COMPUTATIONAL MESHES 

Continuous systems described by partial differential equations respond to differences 
between the state of the system at one spatial location and the state at another spatial location 
that is only infinitesimally distant. The state of the system is defined on a continuum of points 
in the domain. In a finite analog of such a system, the infinity of points in the continuum is 
replaced by a finite number of discrete points, and the infinitesimal distance is replaced by the 
finite distance between neighboring points. For the purposes of computing the differences 
in the state of the system between these points, near neighboring points can be thought of as 
being linked together to form a mesh that covers the domain. The description of the mesh 
consists of a list of the mesh points and their connectivity. The physical relationships between 
the state of the system at one mesh point and that at all others then defines a finite-dimensional 
set of nonlinear algebraic equations that are the exact equations of motion for the finite 
system. The extent to which the dynamics of this finite-dimensional system approximate 
those of the continuum system determines the accuracy and utility of the approximation. 

The computational description of a continuous, time-dependent system, such as a mag- 
netized plasma, has three components: a continuum model of the system that describes 
the evolution of infinitesimally small volume elements for infinitesimally small intervals 
of time; an approximation to the continuum model that describes the evolution of finite 
sized volume elements for infinitesimally small intervals of time; and a description of how 
these finite-sized volume elements evolve over finite time intervals. In this work we have 
chosen resistive magnetohydrodynamics as the continuum model. This will be described in 
Section 3. The finite temporal description will be given in Section 5. Here, and in Section 4, 
we will discuss finite methods of spatial representation. 

2.1. Structured Meshes 

A structured mesh is one in which a predefined logical structure (or order) is assumed 
to exist. For example, in 2D Cartesian coordinates, a structured mesh consists of a product 
of two sets of mesh arrays (the x and y coordinates), with indices i and j, ordered by 
increasing coordinate value. Two indices are required to identify a mesh point: point (/ , j) has 
coordinates x(i),y(j). The mesh is structured logically so that points (/ + 1 , j ) and (/, j + 1) 
are adjacent to point (/, j). This logical structure is assumed to hold for all points in the 
domain and is implicitly used in constructing the finite-dimensional algebraic equations that 
describe the dynamical evolution of the finite system. Structured meshes form the familiar 
quadrilateral grids commonly used in numerical methods. The boundary of the domain 
naturally consists of curves of the form x = const and y = const. (An irregular domain would 
be built up from unions of such meshes.) As neighboring points are logically connected in 
this way, adding and deleting points affects the indexing of all points in the mesh. 

2.2. Unstructured Meshes 

In contrast to a structured mesh, an unstructured mesh is one that has no predefined 
logical structure. An unstructured mesh consists of a set of arbitrarily ordered points. A 
single mesh index suffices to identify a point. Point r, , having coordinates x/ and y t , and 
point i + 1, having coordinates x /+ i and y /+ i, are not necessarily adjacent. 
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Since there is no pre-defined logical structure, the mesh points are not constrained to 
lie along any predetermined curves; they may be arbitrarily distributed in the domain. 
Neighboring points are then connected by line elements to form a mesh of triangles that 
covers the domain. The mesh points r; form the vertices of the triangles, and the connecting 
lines form the triangle edges. The mesh consists of N v vertices, N e edges, and N s triangles, 
with N v < N s <N e . 

With each triangle s we will associate a point r s . This point identifies the location of 
the triangle in the domain. (The definition of r s is not unique. This will be discussed in 
Sections 2.4 and 2.5.) It is also convenient to define the edges of the triangles as directed 
line segments, or vectors 1/j, connecting point i with point j, i.e., 1 = r 7 — r*. Every 

edge e thus has triangle s + L on the left, and triangle s = R on the right. With each edge 
we also associate a unit tangent vector t e = 1 e /l e , and a unit normal vector n e that points 
from the left side to the right side. These mesh elements are sketched in Figure 1 . 

An unstructured mesh is identified and manipulated by means of primary and secondary 
data sets. The primary data set consists of a list of mesh elements. Secondary data sets 
define the connectivity between the primary mesh elements. For example, for 2D meshes 
the spatial representation consists of triangular elements. The primary data set consists of 
a list of cells (triangles), their vertices, and the edges connecting them. Additional data 
sets consist of cross-indexing information that relate the elements of the primary set. For 
example, an edge-indexed array specifies the indices of the cells to the left and right of an 
edge. Other cell-indexed arrays specify the indices of the three vertices and three edges of 
a cell. 


2.3. Primary and Dual Meshes 

Computational meshes, both structured and unstructured, are used not only to describe 
geometric regions, but also to define differential operators. For the latter purpose, it is useful 
to introduce the concept of primary and dual meshes. For a structured mesh, these are often 
referred to as staggered meshes. An example of a two-dimensional structured, staggered 
mesh is shown in Fig. 2. 

In a finite-volume representation, different control volumes are used for different depen- 
dent variables. For example, some dependent variables are defined on at the vertices of the 
primary mesh; the control volume associated with these quantities is a rectangle formed by 
the dual mesh. Other quantities are defined at the vertices of the dual, or staggered, mesh; 
the control volume associated with these quantities is a rectangle formed by the primary 
mesh. Meshes of this type have been used successfully in MHD simulation [16]. 

The concept of primary and dual meshes can be extended to triangular meshes. In this 
case the primary mesh consists of the triangulation of arbitrarily placed points in the plane. 
The mesh points are the vertices of the triangles. The dual mesh consists of polygons that 
surround each vertex. The vertices of the dual polygons can be chosen in several ways. Two 
choices will be discussed below. When taken together, the primary triangular mesh and the 
dual polygon mesh are the generalization of structured, staggered meshes. An example of 
a triangular mesh and its polygon dual are sketched in Fig. 3. 

2.4. The Barycenter, or Centroid, Dual Mesh 

We use a Delaunay triangulation [17] for the primary mesh. The Delaunay triangulation 
maximizes the minimum angle of the triangulation; i.e., of all triangulations of the set P 
the Delaunay triangles are the closest to being equiangular, on average. Because of this 
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I 



FIG. 2. Structured, staggered (dual) meshes. 


property, these triangles are chosen as the primary control volumes. We use a dual mesh 
whose vertices are the centroids, or barycenters, of each triangle. If the coordinates of the 
triangle vertices (the points P) are r v , the coordinates of the vertices of the dual mesh are 
given by 


r s = |(r„i + r v2 + r u3 ), s = 1,2, . . . , N s , 


( 1 ) 


l 



FIG. 3. Triangular (primary) and polygon (dual) meshes. 
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where N s is the number of triangles and the r yz are the three vertices of triangle s. This 
dual mesh has the property that the vertices of the polygons are always interior to their 
corresponding triangles. However, the edges of the dual meshes are not orthogonal, as in 
the Voronoi dual mesh [17]. This complicates the calculation of some dependent variables, 
as will be discussed in Section 4. 


2.5. The Third (Toroidal) Dimension 

One goal of this work is to describe magnetohydrodynamics in toroidal fusion systems. 
The geometry of these systems is symmetric about an axis, and is best described in cylindrical 
(r, 0, z) coordinates. We use the unstructured Delaunay triangular mesh and the barycenter 
dual mesh dual to approximate the geometry in the poloidal (r, z) plane. Since the angular 
(0) coordinate is periodic and since the dominant MHD motions in this geometry are long 
wavelengths in this direction, we have chosen a pseudospectral description using fast Fourier 
transforms (FFTs) for this coordinate. The toroidal mesh is thus structured with a uniform 
mesh spacing A 0 — In/N^, where N# is the number of toroidal mesh points; N# must be 
a power of 2. 

The three-dimensional control volume is sketched in Fig. 4. The elemental volume is 
AV S = r s AcpAa s , where A a s is the planar area of triangle s and r s is the radius of the 
triangle centroid. The Pappus-Guldinus theorem guarantees that this formula is exact. 


2.6. Mesh Refinement 

The use of an unstructured mesh allows for new triangles to be added and old ones 
deleted in a relatively easy manner. Here we use the method of direct dynamic refinement 
[18] (DDR). In this method the mesh is refined or coarsened automatically during the 
execution of the algorithm, according to predetermined criteria. New triangles are added 
to the end of the list, and old triangles are deleted and the list shortened. A new triangle is 
added by introducing a new vertex at the centroid of a triangle to be refined. New edges 



FIG. 4. Three-dimensional control volume. 
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Original Grid 



Refinement 



Refinement and 
One Reconnection 



Second 

Refinement 



FIG. 5. Mesh refinement and edge swapping (from [18]). 


connect this vertex to the three vertices of the original triangle. The original triangle is 
thus divided into three, and two new triangles, three new edges, and one vertex are added 
to the lists. The new edges may need to be swapped between the new vertex and the 
opposing vertices of the three neighboring triangles. The circumcenter test [19] is used 
to determine whether or not edge swapping is required. The new triangulation is thus as 
acute as possible. The addition of a vertex and edge swapping are sketched in Fig. 5. 
Triangle deletion is sketched in Fig. 6. Edges of triangles may also be subdivided. The 
triangle centered densities (mass, momentum, energy, and toroidal vector potential) can 
then be distributed over the new triangles in a conservative manner. Edge centered quantities 
(poloidal vector potential) are interpolated to the new mesh using the methods described in 
Section 4.3. 

Before adaption can occur a triangle must be identified for refinement or coarsening. We 
have used several criteria for this purpose. Different criteria may be useful for different 
problems. One criterion is based on a modified version of the classic interpolation estimate 
originally developed for steady-state hydrodynamic computations [20, 21]. For each triangle 
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Original Grid Point Removal 




Constructing of 
New Cells 


Grid after 
Reconnection and 
Relaxation 


FIG. 6. Mesh coarsening (from [18]). 


s, we compute the normalized second derivative, 

E z 2 s Wu^\ 

s 8 s \VU s \+e\U s \' 


(2a) 


where S s = y/Aa s , U s is any triangle-centered dependent variable, s is a constant between 
0 and 1 (typically, s — 0.2), and the overscore indicates an average over triangle s and its 
three neighbors. (For three-dimensional problems, the maximum of E s over the toroidal 
dimension is taken.) In (2a), the quantity E s is dimensionless and bounded, so that it can be 
used for a variety of problems and dependent variables. Another criterion is based on the 
normalized average gradient: 


&s\VUs\ 

f^max 


(2b) 


For either criterion (2a) or (2b), all triangles for which E s > E R are refined, while all 
triangles for which E s < Ec are coarsened. 

We have found that the proper choices of Er and Ec are as much an art form as they are 
anything else. Choosing too small a value for E R can result in too much refinement and, 
hence, too many triangles. Refinement is done on numerical noise rather than on physical 
structures. Choosing Er too large can result in patchy refinement and poor resolution. The 
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proper choices of Er and Ec also seem to be problem dependent. Much trial and error is 
required. A rule of thumb is to start with Er = 0.9 and Ec =0.1, and make small variations 
until a satisfactory grid results. 

As mesh refinement and coarsening are time-consuming operations they are not done 
every timestep. Mesh refinement is done after every Nr timesteps, and coarsening is done 
after every Nc refinement steps. Values of Nr from 5 to 20, and Nc from 1 to 4, are typical. 
Examples of refinement and coarsening are given in Section 6. 

3. APPROXIMATION TO SPATIAL DIFFERENTIAL OPERATORS 

The differential operators that appear in fluid-like equations are the gradient of a scalar, 
the divergence of both a vector and a tensor, and the curl of vector. We now proceed to define 
approximations to these operators on the triangular, unstructured mesh. We use the method 
of finite volumes as applied to the three-dimensional volume element shown in Fig. 4. 

Consider the triangle in the poloidal (r, z) plane shown in Fig. 7. We define normal and 
tangent unit vectors n e and t e at each edge such that 

= t e X n e , (3) 

where is the toroidal unit vector. (Note that points “into” the page.) The normal and 
tangent unit vectors are given by 

\ e = Ar e e r Az e e z = A/^C (4) 

A z e £ r Ar e e z 


and the area of the triangle is 

xl 2 | = I|l 2 xl 3 | = pl 3 xl 1 |. (6) 


^3 



FIG. 7. Poloidal projection of control volume. 
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The unit normal vector n e points from the left side of edge e (triangle s = L e ) to the right 
side of edge e (triangle s = R e ). 

The finite volume method is used to obtain the approximations to the differential oper- 
ators. In this method differential operators are defined in terms of their integral relations. 
We assume all functions are of the form 

f z, t) = f n (r, z, t )e in ^ (7) 

n 

and then we integrate the appropriate identity over the three-dimensional control volume 
shown in Fig. 4. (Since the toroidal representation is spectral and not finite-difference, the 
limit of the resulting expressions as A0 — > 0 is taken.) This technique assures that the same 
integral relationships are obeyed by the finite difference approximations and their equivalent 
differential operators. 

To obtain an approximation for the gradient of a scalar we substitute Eq. (7) into the 
integral identity 


j VfdV = fnds, 

and use second-order approximations to the volume and surface integrals to obtain 


( 8 ) 


(V/)s = — 1 — y2r e Al e n e f e - — e r + — / s e 0 . 

" ' ■* r r 


r.Aa c 


e=l 


(9) 


The sum is taken over the three edges of triangle s, and the radius of edge e is r e — 
(r ve i + r ve 2 )/ 2 , where r ve \ and r ve2 are the radial coordinates of the vertices connected by 
edge e. The quantity f e is the simple average f e = (f Re + f Le )/ 2, where the values f Re and 
f Le are the values of f s in the triangles lying to the right (R e ) and left (L e ) of edge e. 

Similarly, for the divergence of a vector we use the identity 


V • A dV 


n • A dS 


( 10 ) 


to obtain the approximation 


(V • A) s 


1 

r s Aa s 


3 

E a i a iri 
f eAl e w. e • A e T A(j ) S , 
r s 


( 11 ) 


for the curl of a vector we use 


V x AdS = (hi - Adi 


( 12 ) 


to obtain the approximation 


(V x A) ne 


1 

7~r (r v+ A(f) V + 

r e Al e 


— r y-Affry-.) H A e • t e , 


1 

A a s 


3 

y^A e -i e . 

e=\ 


(13) 


(V x A) 0i = 


(14) 
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Here we have taken surface and line integrals over the faces of the control volume and their 
respective bounding edges. To approximate the divergence of a tensor we use 


V T dV = (b n T dS 


(15) 


to obtain 


1 _ ' 

(V • T), = y— 22r e Al e [e r (n e ■ T e ) r + e 0 ( n e • T e ) 0 + e z (n e • T e ) 2 ] 

5 5 e=\ 

in 1 

H [f*rT(f)rs T ^(})T(f)(f) S + C zT(f)zs\ T [ £(pT(J)rs ^rT<p(f)s\' (lb) 

r s r s 

It is easy to verify from Eqs. (11), (13), and (14) that V • V x A = 0 for these finite 
operators. This is a direct result of the use of consistent integral relations to obtain the finite 
approximations. 


4. THE MHD EQUATIONS: PLACEMENT OF THE VARIABLES ON THE MESH 


In this work we solve the equations of resistive magnetohydrodynamics (MHD). In a 
convenient nondimensional form, they are 



E = — v x B + rjJ/S 
B = V x A 


-V-T 


J = VxB 
dp\ 

T = PVV-BB + \(p + B 2 )l 

g = -V(pv) 


dll 

Yt 


— —V • F 


u — pC -{- 


Y~ 1 


pv 2 H — p ) v + 2 E x B, 

Y ~ 1 


(17) 

(18) 

(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 


where S is the Lundquist number, jj is the resistivity, y is the ratio of specific heats, u is the 
total energy density, F is the energy flux, I is the unit tensor, and T is the Reynolds-Maxwell 
stress tensor. All other quantities have their usual meanings. Note that we have chosen the 
conservation form of the equations. We have found this form satisfactory, even for plasmas 
where the magnetic energy dominates the internal energy. 

Equation (22) is often modified by the addition of an artificial viscosity of the form v Vpv, 
which leads to a vector Laplacian operator on the right-hand side of Eq. (21). This term 
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is necessary to avoid the inevitable buildup of energy at the shortest available wavelength 
(X ~ 8 S = y/ Aa s ) due to nonlinear mode coupling. We have found it necessary to choose 
the viscosity component v such that the cell Reynolds number R A — v8 s /v is of order unity. 


4.1. Hydrodynamic Variables 

The boundary of the poloidal computational region is formed by triangle edges. We define 
the momentum density pv, the energy u and the mass density p at the triangle centroids 
r^. These quantities thus represent the momentum, energy, and mass per unit volume in a 
triangular cell. (The quantities p s AV s , (p\) s A V s , and u s A V s are the total mass, momentum, 
and energy in cell s.) Velocities in a cell are given by \ s = (pv) s / p s - The rate of change of 
these quantities given by applying the differential approximations defined in Section 3.2 to 
Eqs. (17)-(26). For example, the rate of change of mass density in triangle s is given by 


9a 

dt 


1 

r s Aa s 


3 

E in 

r eAl e w e • T F <ps ■> 

. r s 


(27) 


where 


* F e — p e v ne (28) 

is the poloidal mass flux across edge e, and 

F<ps — Ps^cps (29) 

is the mass flux in the toroidal direction. The quantity v ne is the normal component of 
velocity at edge e and is defined as 

Vne = ' (VLe + V Re ) ■ (30) 

Expressions similar to Eqs. (27)-(29) hold for the momentum equation (Eq. (21)), and the 
energy equation (Eq. (24)). 

The advective flux at an edge e is computed using the full donor cell method. For example, 
the right-hand side of Eq. (28) is evaluated as 


( F ne ) a d v — PLe^nei if Vne > 0, (31a) 

(Fne)cidv — PRe^ne-> if Vne < 0? (21b) 

where L e and R e are the indices of the triangles to the left and right of edge e 9 respectively. 
This method introduces numerical diffusion of order v ne 8 s / 2, where 8 S ^ Aa s . While 
this technique may be too diffusive for highly accurate shock calculations, it is quite adequate 
to describe the relatively slow motions of interest in tokamak dynamics. Problems involving 
strong shocks may require a higher order treatment. 

4.2. Electromagnetic Variables 

The primary electromagnetic variable in this formulation is the vector potential A. We 
define A r and A z at the triangles edges e, and A $ at the triangle centroids s. Then Eqs. (13) 
and (14) define B ne , the component of B in the poloidal plane normal to a triangle edge, 
and the toroidal component of B at the triangle centroid. (Note that (V • B ) s = 0.) 
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Polygon edge p e 
Triangle edge e 


FIG. 8. Triangle and polygon edges. 


The procedure described above defines only the component of B normal to each triangle 
edge. To uniquely determine the magnetic field we must also define another independent 
component of B in the poloidal plane. This is done by integrating Eq. (12) over the surface 
of the dual polygon p e that crosses an edge e , as shown in Fig. 8. The polygon edge has a 
unit normal vector n pe and a unit tangent vector t pe . This defines B npe , the components of 
B normal to the polygon edge. The cylindrical components of the poloidal field at a triangle 
edge are then given by 


1 

Bye — ~^\B ne Ti pez B npe n ez ) (33) 

1 

B ze — ~^(yB npe n er B ne n per ), (34) 


where 


A = e 0 • (n e x n pe ) / 0, 


(35) 


from which the tangential component of B at edge e is computed as 

Bte — B ze n er B re n ez . (36) 

Similar relationships hold for the current density J. (Note that if the mesh consists of 
Delaunay triangles and Voronoi polygons the dual meshes are orthogonal and this calculation 
is simplified.) 

In light of Eq. (17), we define the electric field E at the same spatial locations as the 
vector potential A. The normal and tangential components of the electric field at a triangle 
edge are given by 


Ene — V(j)eBte T V neBfe T T)J ne /S, (37) 

Efe — V ne B(f)e T VcfieBfje YjJfy/S. (38) 

The toroidal electric field at the triangle centroids is given by 


E(ps — VzsBys T Vf$ B zs T tj J (j) S / S . 


(39) 
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We note that Eqs. (11) and (16) are the discrete equivalent of Gauss’ theorem, and 
Eqs. (13) and (14) are the discrete equivalent of Stokes’ theorem. Thus, all conservation laws 
obeyed by the resistive MHD equations have equivalent discrete statements. For example, 
the rate of change of toroidal flux through the surface of any triangle is given by the line 
integral of the tangential electric field around the edges of the triangle, in exact equivalence 
to Lenz’s law. This is a direct result of Eqs. (17) and (19) and the definition of the differential 
operators given here. It is independent of the form of Ohm’s law. If the surface is constructed 
from several triangles, the line integrals over the interior edges cancel, and only the line 
integral of the tangential electric field around the bounding surface remains. If that surface 
is a perfect electrical conductor, the tangential electric field vanishes on the boundary and 
the total toroidal flux is conserved, even in a resistive plasma. 

4.3. Averages and Interpolation 

In Eqs. (37)-(39), an overscore indicates that an average should be taken or that inter- 
polation be performed. Several types of interpolation are required in the present algorithm, 
especially during mesh refinement and coarsening. These are discussed in this section. 

Interpolation from triangle centroids to edges is a simple average between adjacent tri- 
angles: 


fe = \ifRe + he)- (40) 

Interpolation from vertices or edges to triangle centroids is also a simple average. For 
functions defined on vertices, 


fs — |(/u 1 + fv 2 + fv 3), (41) 

and for functions defined at edges, 

fs = life 1 + fe 2 + fel), (42) 

where vl, v2, v3, and el, e2, e3 are the three vertices and edges of triangle 5 , respectively. 
Interpolation from edges to vertices is given by 

( 43a > 

ve e> 

where N ve is the number of edges meeting at vertex v, and e indicates a sum over those 
edges, while the interpolation formula from vertices to edges is 

fe = \ifv\+ fvl), (43b) 

where vl and v2 are the two vertices joined by edge e. 

For interpolation from triangle centroids to vertices, we use a pseudo-Laplacian weighted 
average [22]. In this approach, the interpolated value of a function at vertex v is given by 
the weighted average 


7 _ E,<(i + mv )jf 

Tv y,A i + «v) ’ 


(44) 
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where the prime (s') indicates that the sums are taken over all triangles sharing vertex v, 
and not over all triangles N s . 

The weights w s are determined by requiring that they be as small as possible and that the 
interpolation be exact for linear functions. We then minimize the functional 

F(w s ) = Y j u>1 ( 45 ) 

s' 

subject to the constraints 

L, = + w s')(r S ' — r v ) = 0, (46a) 

s' 

L z = y](l + uv)0v - Z„) = 0, (47b) 

s' 

where (r v , z v ) and (r s , z s ) are the coordinates of the vertex and the centroids. The result is 
that 

w s = X r Ar s + X Z A z s , (48) 

where Ar s = r s — r v , A z s = z s — z v , and X r and X z are Lagrange multipliers given by 


^ _ Rzhz — Rrhz 
r ~ II — 1 2 ’ 

l rr l zz l rz 

(49) 

R z I rz Rzhr 
kl ~ II -I 1 ’ 

L rr l zz l rz 

(50) 

< 

II 

(51) 

*«o 

N 

<1 

II 

(52) 

hr = yyAr s 0 2 , 

s' 

(53) 

izz = y>^) 2 , 

(54) 

lrz = ^ Ar s ' Az s ' . 

(55) 


s' 


Equations (40)-(44) have low order of accuracy. The use of higher order interpolation 
methods, especially in place of Eq. (40), can be shown to lead to a non-Hermitian formu- 
lation of sound waves and resulting unphysical behavior. 

A complication is that neither Eqs. (40) and (42), Eqs. (41) and (44), nor Eqs. (43a) and 
(43b) are exact inverses of each other. Thus, for example, interpolation from centroids to 
vertices using Eq. (44), followed directly by interpolation from vertices to centroids using 
Eq. (41), introduces errors. Heuristically, these errors do not seem critical to the results 
obtained with the algorithm, but their affect on problems in parameter regimes other than 
those studied in the present work cannot be assessed. 
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4.4. Boundary Conditions 

Since the computational boundary consists of triangle edges, the mass density in triangles 
adjacent to the boundary is completely determined by the normal velocity v ne at the bound- 
ary. The momentum density and energy density also require that the normal component 
of the magnetic field, B ne , and the total pressure, p + B 2 , be specified. For a nonporous, 
perfectly conducting wall, the appropriate boundary conditions are that B ne — v ne = 0, and 
that the normal gradient of the total pressure vanish. Implementation of boundary conditions 
is aided by introducing ghost triangles that lie outside the boundary and are reflections of 
interior triangles that contain a boundary edge. For the electromagnetic variables it is suffi- 
cient to specify the electric field tangent to the boundary. Thus, for a perfectly conducting 
wall E te — Efa = 0, where E$ e is the average of the toroidal electric field in a boundary cell 
and its reflected ghost cell. 

When the inner boundary extends to r = 0, regularity conditions must be applied to 
the complex Fourier coefficients. For scalars ( pv z , A z , B z , J z , p or w, and p), all Fourier 
components vanish except n = 0, which have zero derivative. For vector components 
(pv r , pV(j), A r , A#, B r , B^, J, , and J#), all Fourier components vanish, except n — 1, which 
have zero derivative. However, for n — 1 this condition is only directly applied to the 
r-component of vectors. The ^-component is found from the relation V r + / V# = 0, which 
assures that the Cartesian (x,y) components of the vector are unique at r — 0. 

Boundary conditions corresponding to applied tangential or toroidal electric fields and 
to supersonic inflow (specified upstream velocity, pressure and density) and outflow (zero 
normal derivative of pressure, density, and velocity) have also been implemented. 


5. TIME INTEGRATION 

As is appropriate for sound and Alfven waves the time integration algorithm uses an 
explicit leapfrog method with predictor-corrector steps to stabilize the nonlinear advective 
terms. The leapfrog scheme is second-order accurate for uniform time steps, while the 
predictor-corrector is formally first-order accurate in time. The velocity and momentum 
are defined at time t n . The energy density, mass density, and vector potential are defined 
at time C +1/2 . The time step can be arbitrarily large; the semi-implicit method [13-16] is 
used to remove the CFL time-step restriction and is second-order accurate for uniform time 
steps. Artificial viscosity, which is required for nonlinear numerical stability, is treated fully 
implicitly and is, therefore, first-order accurate. 

The time advance proceeds by means of operator splitting, i.e., 


9U 

au 

au 

au 



— — 

+ — — 

— 


dt 

total dt 

explicit 

O f 

semi-implicit 

viscous 
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(57a) 
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where U is the state- vector describing the system, and F exp ii cit , F sem i-im P iicit, and F v i SC ous 
represent the explicit (wave-like and advective), semi-implicit, and viscous terms that appear 
on the right-hand side of the equations. Details of these methods are given in the following 
sections. 


5.1. Explicit Advance 

Wave-like and advective terms are advanced explicitly with At chosen for accuracy and 
computational convenience rather than numerical stability. The explicit part of the algorithm 
is 


pv* — p\ n 1/2 
At 


— V • (pvv)" 1/2 


p\ n + 1 / 2 — py n ! / 2 
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(r - 1 )[«* - ( pv 2 ) n+I/ 2 - ( B 2 r +i ] 

-x n+ 1/2 x B" +1 + r]J n+l /S 
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(58) 

(59) 

(60) 
(61) 
(62) 

(63) 

(64) 

(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 

(70) 


Total mass, momentum, and magnetic flux are exactly conserved. Because the pressure, 
magnetic field, mass density, and momentum are defined at different time levels, the sum 
of the kinetic, magnetic, and internal energies is exactly conserved in the limit At -> 0, 
independently of spatial discretization. (The volume integral of the quantity u is exactly 
conserved independently of At.) The predictor-corrector steps introduce an additional 
diffusion of order vAt/2 that can exceed the diffusion from the donor cell fluxes when the 
time step exceeds the explicit CFL stability limit. 
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5.2. Semi-Implicit and Implicit Solutions 

We use the semi-implicit method [13-16] to remove the CFL time step restriction for 
numerical stability associated with the explicit advance described in Section 5.1. This 

restriction is of the form C At/8 <1, where C = \Jv 2 + Cj + V\ is the characteristic 
speed for the propagation of disturbances, and <5 is a measure of the linear size of a zone 
(here proportional to the square root of the triangle area). With the use of the semi-implicit 
method the algorithm becomes numerically stable at arbitrary At so that the time step can be 
chosen for reasons of accuracy or computational convenience rather than numerical stability. 
The time step remains limited by the advective CFL stability condition VAt/8 < 1, where 

V is the local flow speed. This is not a significant restriction when V / y C] + V\ 1, as is 

the case for many fusion applications. When V/y C] + V\ ^ 1, as is the case for shocks, 
the algorithm becomes explicit. This restriction can thus be viewed as an accuracy condition. 

In this work we use a simple vector Laplacian semi-implicit operator [16]. This term is 
added to and subtracted from the right-hand side of the momentum equation at the new and 
old time levels. The semi-implicit advance is 

(1 - a A tV 2 )(pv)** = Gov)* - aAtV 2 (p\) n , (71) 


where a is the semi-implicit coefficient given by 

~ ' At \ 2 


ki At 


\ AfcFL , 


- 1 


for At > A^cfl> 


a = 0 for At < Af C FL, 


(72) 

(73) 


where Gov)* is the value of momentum obtained from the explicit advance (Eq. (59)), 
& max ~ l/<$y is the largest wave-number resolved on the mesh, Ar C FL is the maximum time 
step allowed by the CFL restriction for normal modes, and a is a constant > 1 . Since Eq. (7 1 ) 
is second order in space, boundary conditions must be applied to the tangential (</> and z) 
components of the momentum density. For inviscid flow, the normal derivative is set to zero 
at solid boundaries; for viscous flow, the tangential momentum is set to zero. In either case, 
the normal component vanishes at solid boundaries. For supersonic inflow, the upstream 
velocity is specified. For outflow, the normal derivatives of all components are set to zero. 
Boundary conditions at r — 0 are as described in Section 4.4. 

The time step is completed with the implicit viscous advance 

(1 - vA?V 2 )(pv)” +1 = Gov)**, (74) 


where v is a (possibly spatially dependent) artificial viscosity coefficient. All components of 
momentum are set to zero at solid boundaries. For inflow, outflow, and r = 0, the boundary 
conditions are as described for Eq. (71). 

The vector Laplacian operator appearing in Eqs. (71) and (74) requires the definition of 
the scalar Laplacian. This is accomplished by the successive application of the gradient and 
divergence operators defined in Eqs. (9) and (11). When combined with the condition 


(V/) G = (V/)*, 


(75) 
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FIG. 9. Matrix structure for the Laplacian operator on an unstructured mesh. 


where the subscripts G and B represent values in ghost and boundary triangles, respectively, 
the resulting operator is self-adjoint. 

Since the mesh is unstructured, the N s x N s matrices corresponding to the operators 
appearing in Eqs. (71) and (74) are not banded but are sparse. An example of the structure 
pattern for a case with 320 triangles is shown in Fig. 9. Matrix inversion is performed with 
a conjugate gradient (CG) algorithm with diagonal preconditioning [23]. Since this method 
is iterative, the full N s x N s matrix never needs to be stored. Good convergence properties 
have been found, even with relatively large time steps. 

At this time the resistivity is treated explicity. Since ?]/5 < 1 we have not found this to 
be computationally restrictive. 


6. APPLICATIONS 

The algorithm described above has been applied to several nonlinear test problems, both 
two- and three-dimensional. The code based on the algorithm is called TRIM, for TRIangular 
MHD. The application of TRIM to these test problems is described in the following sections. 

6.1. The Hydrodynamic Shock Tube Problem 

A standard problem for testing hydrodynamic algorithms has been defined by Sod [24]. 
The initial conditions consist of two fluids with different uniform properties separated by 
a membrane. The fluid to the left of the membrane has both pressure p L and density pi 
equal to 1 . The fluid to the right of the membrane has pr =0.1 and pR = 0. 125. The initial 
velocity is zero and the ratio of specific heats is y = 1.4 (air). The magnetic field is zero. 
These conditions are sketched in Fig. 10. 

At t = 0 the membrane is ruptured and the fluid reacts dynamically. This Riemann 
problem is one of the few fully nonlinear problems that has a known analytic solution [9], 
and is therefore valuable for testing numerical algorithms. The solution consists of an 
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0 L 0.5 R 1 


PL - 1 PR - ai 
PL = 1 PR = °- 125 

FIG. 10. Initial conditions for hydrodynamic shock tube problem. 


expansion wave traveling to the left and a shock wave and a contact discontinuity traveling 
to the right, all with known velocities. 

We have applied the TRIM algorithm to this problem. The time integration is explicit 
and the artificial viscosity v is set to zero. While this test problem is one-dimensional, the 
triangular grid in TRIM requires that a two-dimensional problem be solved. The mesh is 
shown in Fig. 1 1 . It contains 1250 triangles. In this figure, the initial membrane is horizontal 
and centered at z = 0.5. As the solution proceeds in time no spatial variation develops in the 
direction parallel to the membrane. The solution thus remains one-dimensional, even with 
the two-dimensional algorithm. The analytic solution at t = 0.1 is shown in Figs. 12a-c. 

The results of TRIM with the mesh shown in Fig. 11 is shown in Figs. 13a-c at 
t = 0.1. The magnitude of the pressure and velocity in the region between the shock and 
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1001.0 


FIG. 11. Mesh for hydrodynamic shock tube problem, with superimposed contours. 
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FIG. 12. (a-c) Analytic solution to the hydrodynamic shock problem. 


the expansion fan are quite accurate. (Note the because of the normalization the pressure 
in TRIM appears to be twice the pressure in the analytic solution.) As is anticipated, the 
numerical diffusion introduced by the first-order upwind treatment of the interface fluxes 
has resulted in a considerable smoothing of the discontinuities. This is especially notice- 
able in the density. The contact discontinuity, which is an interface separating regions of 
different density but equal pressure and velocity, has been considerably smeared out. This 
structure is particularly difficult to treat numerically. In contrast with a shock, there are 
no nonlinear processes that continue to generate a contact discontinuity in opposition to 
numerical diffusion; it is merely an interface between two states of different density. The 
effect of any diffusion in the algorithm is felt most strongly here. 

One solution to the problem of numerical diffusion is to employ a higher order approxi- 
mation to the interface fluxes. Another solution is to use a low order method but to reduce 
the diffusion by adaptively refining the mesh in the regions near the discontinuities. We 
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FIG. 13. (a-f ) Comparison of numerical solutions without (a-c) and with (d-f ) dynamic mesh adaption. 

have applied the mesh refinement techniques described in Section 2.7 to this problem. 
For this problem we used the mass density p as the dependent variable in the refinement 
criterion, Eq. (2a), with s = 0.2, E R = 0.8, and Eq = 0.2. (Smaller values of E R resulted 
in an unacceptable number of triangles; larger values resulted in patchy refinement. See 
Section 2.6.) Both refinement and coarsening were done every 50 timesteps. The initial 
mesh is shown in Fig. 14a. This mesh has been refined in order to initially capture the 
pressure and density discontinuities. The adaptively refined mesh at t = 0.1 is shown in 
Fig. 14b. The algorithm has adapted the mesh to the dynamically evolving shock, contact 
discontinuity, and expansion front. The initial mesh had 7777 triangles, and the dynamically 
evolving mesh contained up to 34,415 triangles. 

In Figs. 13a-f we compare the solution at f = 0.1, with and without adaption. All 
features are sharper with dynamic mesh refinement than without. In both cases, the jump 
conditions are exact. The discontinuities become sharper as the grid is refined. However, we 
are unable to furnish an estimate of the error as a function of grid refinement, as this is the only 
satisfactory grid we could produce. As anticipated, the major error is in the contact discon- 
tinuity. This error can only be eliminated by the use of a higher order method, such as a 
Riemann solver. This is beyond the scope of the present work. We emphasize that problems 
involving strong shocks are uncommon in fusion plasmas, so that low-order methods are 
sufficient for these applications. 

6.2. The Magnetohydrodynamic Shock Tube Problem 

The hydrodynamic shock tube solution described in the previous section has been ex- 
tended to MHD by Brio and Wu [25]. The thermodynamic properties of the left and right 
states are the same in the purely hydrodynamic case. A uniform magnetic field B x is imposed 
in the direction (x) perpendicular to the membrane. The component of the magnetic field B y 
parallel to the membrane is discontinuous at the membrane, with B yL = 1 and B yR = — 1 . 
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FIG. 14 . (a) Initial mesh for dynamic adaption; (b) dynamically adapted mesh at t = 0.1. 
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FIG. 15. Initial magnetic field vectors for the MHD shock tube problem. 


The membrane is thus a current sheet in the z-direction. The magnetic configuration is 
sketched in Fig. 15. 

The dynamics after the membrane is ruptured are much more complex than in the purely 
hydrodynamic case; we refer the reader to Ref. [25] for details. In Fig. 16 we present 
our two-dimensional solution of this problem. This can be compared with the more finely 
resolved one-dimensional solution of Brio and Wu [25]. We find that most of the details 
of the Brio-Wu solution are reproduced in our results, although the effect of the low-order 
diffusion is again apparent, especially near the contact discontinuity. We have also repeated 
the calculation with the component of magnetic field parallel to the membrane rotated by 
7r/2 and find identical results for this polarization. 

We have also applied mesh refinement and coarsening to the MHD shock tube prob- 
lem. The refinement and coarsening criteria are the same as described in Section 6.1. Fig- 
ure 17 shows a comparison of the mass density p with and without dynamic mesh adaption. 
Finer structure is observed when dynamic mesh refinement is implemented. However, the 
contact discontinuity is still poorly resolved due to the low order calculation of the interface 
fluxes. 


6.3. Supersonic Flow Past a Sphere 

We now consider the case of steady axisymmetric supersonic flow past a solid sphere. 
The sphere is centered at r = z = 0. The boundary of the sphere is at r 2 + z 2 = 1. 
The outer boundary of the computational domain is a cylinder of radius r = 10 with 
ends at z = ±10. Supersonic inflow conditions are imposed at the upper boundary (z = 10) 
and outflow conditions are imposed at the lower (z = — 10) and outer (r = 10) boundaries. 
Regularity conditions are imposed at r — 0. 

The initial conditions consist of uniform pressure, density, and axial velocity (v z ). The 
initial radial velocity (v r ) is set to zero. The parameters are chosen so that the initial axial 
velocity corresponds to an upstream Mach number M = v/c s = 2. 

With the magnetic field set to zero, the hydrodynamic equations are integrated forward 
in time until a steady (time independent) state is reached. A small amount of viscosity 
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Pressure 







FIG. 16. Numerical solution of the MHD shock tube problem. This can be compared with the solution given 
in [25]. 


No Dynamic 




FIG. 17. Comparison of the results for the MHD shock tube problem without and with dynamic mesh adaption. 
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FIG. 18. Velocity vectors and pressure contours for M = 2 flow past a sphere. 


(v = 10 -3 ) was used in this calculation. In the steady state a shock forms in front of the 
sphere as the flow is deviated around it. The flow separates from the surface of the sphere 
and forms a recirculating eddy behind it. The pressure contours and velocity vectors in the 
vicinity of the sphere are shown in Fig. 18. 

This calculation was performed using an adaptive mesh along with the criterion given by 
Eq. (2b). The mass density was used to compute the normalized gradient. The mesh was 
locally refined if E s > 0.95, and coarsened if E s < 0.05. Refinement and coarsening were 
performed every 500 timesteps. The initial and final meshes are shown in Figs. 19a and b. 
The final mesh has 13,360 triangles. 

We have repeated this calculation, but with a dipole magnetic field embedded in the 
sphere at t = 0. The axis of the dipole is aligned with the upstream flow direction. This 
field is given by the vector potential 


A(t> B ° ( r 2 + z 2 ) V 2 ’ 


(76) 


where Bo is the strength of the magnetic field at r = 1 , z = 0. The resistivity was taken 
to be uniform and to correspond to a Lundquist number of S = 10 3 . A coarse mesh was 
used. The mesh in the vicinity of the sphere is shown in Fig. 20. The initial magnetic field 
is shown in Fig. 21. 

As the calculation proceeds, the initial dipole field remains embedded in the spherical 
boundary and is swept downstream behind the sphere. The shock remains in front of the 
sphere. The final magnetic field for the case Bq = 0.5 is shown Fig. 22, and the adapted 
grid is shown in Fig. 23. 
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FIG. 19. (a) Initial mesh; (b) final, adapted mesh for M = 2 flow past a sphere. 


6.4. Nonlinear Evolution of Toroidal Instabilities 

One of the standard applications of MHD to fusion plasmas is the linear growth and 
nonlinear saturation of instabilities. These instabilities can occur because toroidal equilibria 
are not necessarily minimum energy states, even though they are extrema of the energy. 
Equilibria that are local maxima of the energy are unstable, with small deviations from 
the initial state growing exponentially in time. Determining the stability of equilibria is an 
important problem in the design of a fusion experiment. Even stable equilibria can be driven 
unstable by diffusive processes [7]. 

It is to be emphasized that algorithms of the type described in this paper are not the most 
efficient or accurate way of computing linear stability. Specialized algorithms that find the 
eigenvalues of the linearized MHD operator are better suited to this problem [26]. Nonethe- 
less, linear stability problems are among the few three-dimensional toroidal problems with 
known solutions (generally obtained with the specialized algorithms [26]) and are therefore 
valuable benchmarks for nonlinear, time-dependent algorithms. 

Physically, linearly unstable MHD modes are of interest only if they impart some obser- 
vable, and hence, finite and nonlinear, perturbation to the physical system. The details of 
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FIG. 20. Mesh used for flow past a sphere with an embedded dipole magnetic field. 


the nonlinear evolution of linearly unstable modes requires that algorithms of the present 
type be employed. 

6.4. 1 . Solov’ ev Equilibrium 

Axially symmetric force balance in a torus is given by solutions to the Grad-Shafranov 
equation 




9 Vt 1/ 1 ? dP dF 

r 2 V • — = -r 2 F , 

r 2 2 di]f dxfs 


(77) 


where z) = r A ^ is the poloidal flux, and the pressure P (t/ 0 and the toroidal flux 
function F (ifr) — rB ^ are arbitrary functions of x/f. An analytic solution has been given by 
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FIG. 21. Initial dipole magnetic field. 


Solov’ev [27], With 


4 1 +K 2 

W)-l — 
e l k 1 

2b 

= 1-VO 

£K 


the poloidal flux and toroidal field are 

f(r, z) = \ 
£ l 

c 


(r 2 + b 2 )z 2 


= - + ©(*), 


(f- 1) 

(78) 

/2 + c 

(79) 

(r 2 - l) 2 " 


+ 4 J 

(80) 


(81) 


where e = a/R is the inverse aspect ratio, k is the elongation, b is a diamagnetic factor, 
and C is a normalization constant that determines the strength of the vacuum toroidal 
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FIG. 22. Final magnetic field configuration for M = 2 flow past a sphere with an embedded initial dipole 
field. 


field. Contours of x/r and P with k = 1 and b — 0 are shown in Fig. 24. Since b — 0, this 
equilibrium has no poloidal current (/, = J z = 0). 

Linear stability and comparison with previous results. The linear stability of the Solov’ev 
equilibrium to large-scale MHD modes has been extensively studied using specialized 
eigenvalue techniques [28]. We have developed a “linearized” modification of the TRIM 
algorithm to study the linear stability of this equilibrium. This modification is possible 
because of the pseudospectral representation. In these calculations the initial conditions 
consist of the axisymmetric ( n = 0) equilibrium, such as Eqs. (80) and (81), along with very 
low amplitude random noise in the velocity field of a single n = no >0 toroidal mode. All 
other toroidal modes are set to zero initially. The calculation then proceeds as described 
in Section 5, except that after each time step the amplitudes of all modes with n / no are 
reset to zero and the n = 0 (equilibrium) component is restored to its initial value. This 
effectively disables any nonlinear or quasilinear mode couplings and affords a good ap- 
proximation to the solution of the linearized equations. The magnetic and kinetic energies 
of an unstable mode will grow exponentially with time as exp(2 yt), allowing the growth 
rate (eigenvalue) y to be calculated. The self-similar spatial structure of the growing mode 
defines the eigenvector. 
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FIG. 23. Adapted grid for flow past a sphere with embedded dipole magnetic field. 


The Solov’ev equilibrium can be completely parameterized in terms of the three non- 
dimensional constants s (the inverse aspect ratio), k (the elongation), and qo (the safety 
factor on axis). In Ref. [28] values of normalized growth rate yiA (where Ta is the Alfven 
transit time) were obtained over a range of qo for values of s = 1 /3 and k = 1 and 2. Here 
we have primarily focused our attention on the cases with qo = 0.5, which exhibit robust 
instability to ideal MHD modes for these values of s and k . Special cases with qo — 0 and 
qo — 0.8 will also be described. Also, we have used a boundary condition that corresponds 
to a perfectly conducting boundary placed at the plasma edge. Thus, only internal (rigid 
boundary) modes are considered. 
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Poloidal Flux Pressure 




FIG. 24. Contours of poloidal flux and pressure for the Solov’ev toroidal equilibrium with k — 1. 


When the semi-implicit method is used to take time steps that exceed the CFL limit, the 
growth rate will be a function of the time step [16]. This is discussed in more detail later in 
this section. The actual growth rate is obtained in the limit At — > 0. (In practice, we only 
take the limit At/AtcFL — ► 1.) Here we have performed such a convergence study only for 
the two cases e — 1 /3, no = 2, qo = 0.5, and k = 1 and 2. Quantitative comparisons should 
be made for these cases only. 

For the converged cases (qo = 0.5, no = 2), we find for k — 1 , y x A = 0. 122, as compared 
with yx A = 0. 158 from Ref. [28] (when converted to our normalization); for k — 2, we find 
yx A =0.22, compared with yx A =0.28 from Ref. [28]. The time-step converged growth 
rates determined from TRIM are about 20% lower than those of Ref. [28] and are consistently 
lower for the nonconverged cases. A quantitative result for the case with qo — 0, k = 1, is 
not given in Ref. [28], where this mode is identified as “an m — 0 mode and not a kink.” 

The eigenfunctions (spatial structure) for the cases no = 2, 3, and 4, qo = 0.5, and k = 1 
are shown in Figs. 25-21. In each figure, velocity vectors of the real part of the poloidal 
velocity (v r , v z ) and contours of the imaginary part of the toroidal velocity (i^) for mode 
no are shown. For the no — 2 mode, the poloidal structure is dominantly m — 1, while for 
the no = 4 mode, the poloidal structure is dominantly m = 2. This is consistent with the 
value qo = 0.5, and in agreement with the results of Ref. [28]. The eigenfunction for the 
case no = 2, qo = 0, k — 1 is shown in Fig. 28. It is easy to see the dominant m — 0 inter- 
change structure of this mode. 

The equilibrium for the case k = 2, qo = 0.5 is shown in Fig. 29. Velocity eigenfunctions 
for the cases k = 2, qo = 0.5, and no = 2 and 3 are shown in Figs. 30 and 3 1 . The dominant 
poloidal mode structure is in agreement with that of Ref. [28]. Finally, in Fig. 32 we display 
the velocity eigenfunction for the case k — 2, qo — 0.8, no — 1 . The rigid m = 1 displacement 
of the mode is in contrast with the vortex structure displayed by the other unstable modes 
found here. 

Effect of time step on linear growth rate. As discussed in Section 5, TRIM uses the semi- 
implicit method [13-16] to achieve time steps in excess of that set by numerical stability 
constraints. As shown in Ref. [16] the semi-implicit method reduces the characteristic 
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FIG. 25. Linear eigenfunction for the n = 2 mode in the Solov’ev equilibrium with k = l, q 0 = 0.5, Rfa = 3. 
Velocity vectors display the real part of the poloidal velocity (tv, tv), contours display the imaginary part of the 
toroidal velocity (v f ). 


frequency of a mode with wave number k by a factor 



FIG. 26. Linear eigenfunction for the n = 3 mode in the Solov’ev equilibrium with k = 1 , q 0 = 0.5, R/a = 3. 
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FIG. 27. Linear eigenfunction for the n = 4 mode in the Solov’ev equilibrium with k = 1 , q {) = 0.5, Rfa = 3. 


calculation with At = A^cfl- Using Eq. (72), this expression can be rewritten as 



(83) 


where a is a constant of order unity and k mdx is the maximum value of k that can be resolved 
on the grid. For unstable modes, the frequency co becomes the growth rate y. 

For the case no = 2, k =1, and (70 = 0.5 (see Fig. 25), we have k mSLX = 164 (corresponding 
to 2090 triangles), k ^ 15 (corresponding to X r & X z 2a ^ 0.6), and a = 1 .5. This case 
was run with a small amount of artificial viscosity, v = 10 -4 . In Fig. 33 we plot the ratio 
y/Xo obtained from both TRIM and from Eq. (83) as a function of Af/Af C FL for this case. 
These results substantially confirm the effect of the semi-implicit method on the growth 
rate. The discrepancy between the two curves may be due to the artificial viscosity or other 
numerical damping inherent in the algorithm. 

Nonlinear results. The linear results described above were obtained by removing all 
nonlinear interactions, freezing the n = 0 component of the solution, and allowing only a 
single mode with n — no to evolve. To obtain the full nonlinear dynamical evolution of 
the unstable equilibrium, all modes with n > 0 are initially perturbed, all nonlinear inter- 
actions are restored, and the n = 0 component is allowed to evolve under the influence 
of unbalanced forces, resistivity, viscosity, and nonlinear effects. Since these cases may in- 
volve large amplitude displacements and considerable dynamics, larger values of resistivity 
and artificial viscosity are used than in the strictly linear results. Typical values are S = 1 0 4 
and v — 10 -2 . Total toroidal flux is conserved. 

In Fig. 34 we plot the time evolution of the kinetic energy in the modes 1 < n < 5 (corre- 
sponding to 16 toroidal mesh points) for the case k = 1 , qo = 0.5, R/a = 3. The time step 





FIG. 28. Linear eigenfunction for the n = 2 mode in the Solov’ev equilibrium with k = \, = 0, R/a = 2. 

Note the m = 0 interchange structure. 


is such that At / A^cfl = 1 1 - 5. Because of the finite resistivity, the total toroidal current de- 
cays and its profile peaks during the evolution, thus altering the linear stability properties of 
the discharge. The n = 2 mode is linearly unstable, grows to finite amplitude, and saturates. 
The n = 4 mode, which was found to be linearly unstable by both TRIM and in Ref. [28], 
exhibits initial exponential growth at approximately twice the rate of the n = 2 mode. This 


Poloidal Flux 


Pressure 




FIG. 29. Poloidal flux surfaces and pressure contours for the elongated Solov’ev equilibrium k — 2, 
q 0 = 0.5, R/a = 3. 
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FIG. 30. Linear eigenfunction for the n = 2 mode in the Solov’ev equilibrium with k = 2, q 0 = 0.5, R/a = 3. 


and the delayed onset of the mode imply that this mode is driven nonlinearly by the n = 2 
mode rather than by inherent linear instability. The change in the linear stability properties 
of this mode may be due to the modification of the n = 0 component by resistive diffusion, 
but we have not verified this conjecture. The n — 1 and n — 5 modes exhibit complete lin- 
ear stability. The small amplitude increases at late time in these modes is due to nonlinear 
spectral broadening from the saturation of the n — 2 mode. The n — 3 mode shows some 
small indication of linear instability late in the calculation but prior to n — 2 saturation. We 
cannot rule out nonlinear spectral broadening as the cause of this increase. In Fig. 35 we plot 
contours of the pressure in the saturated state at four different toroidal locations spanning 
one-half of the torus. The predominantly n = 2, m = 1 helical displacement of the plasma 
column is evident. 

The same calculation has been performed for the case k — 2, R/a = 3, q 0 = 0.5, with 32 
toroidal mesh points (corresponding to toroidal mode numbers 0 < n < 10 after dealiasing). 
The magnetic energy in the modes 1 < n < 10 is shown in Fig. 36. The modes n — 2 and 
n = 3 exhibit robust linear instability. The n = 4 mode also indicates initial linear instability, 
but makes a transition to being nonlinearly driven by the n = 2 mode later in the calculation. 
The n — 1 and n = 5 modes are driven by the nonlinear interaction of the n = 2 and n — 3 
modes, and the n = 6 mode is nonlinearly driven by the n = 3 mode. All other modes 
appear to be driven by nonlinear spectral broadening. The finite amplitude of all modes 
at the end of the calculation probably indicates that more toroidal mesh points (modes) 
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FIG.31. Linear eigenfunction for the n = 3 mode in the Solov’ev equilibrium with k = 2, q Q = 0.5, Rfa — 3. 

are required for proper resolution of the nonlinear state. In Fig. 37 we plot contours of 
the pressure in the saturated state at four different toroidal locations spanning one-half the 
torus. Again, the helical displacement is topologically m = 1, n = 2 dominant. 

6.4.2. ITER Equilibrium 

In Figs. 38a and b we display poloidal flux and pressure contours for an equilibrium 
that is representative of ITER, an international fusion test reactor that is presently being 
designed [29]. The outer boundary is the separatrix, or last closed flux surface; it is intended 
that the plasma is confined within this surface. This equilibrium been constructed to be 
robustly unstable to an internal kink mode with toroidal mode number n — 2. The safety 
factor profile is shown in Fig. 39. (This equilibrium is unphysical, but it retains the ITER 
cross-sectional shape and is robustly unstable. It thus serves as a demonstration problem.) 
The linear instability has been computed with the GATO code [26], which directly solves 
the resulting linear eigenvalue problem. 

The boundary of the calculation is the separatrix, which we take to be rigid and per- 
fectly conducting. The unstructured mesh inside the separatrix with N s = 5728 is shown in 
Fig. 40. The equilibrium is initialized to this mesh by cubic spline interpolation, and the 
resulting force imbalance is resolved with viscous damping. We have also found it useful to 
introduce spatially dependent resistivity, with S = 10 6 near the magnetic axis and S = 10 4 
near the separatrix. Thus resistive flow is always present and true static equilibrium is not 
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FIG. 34. Kinetic energy in the modes 1 <ri <5 versus time for the nonlinear evolution of the Solov’ev 
equilibrium with k = 1, q Q = 0.5, R/a = 3. 


achieved. The resistivity also causes the current to peak near the magnetic axis, thus altering 
the safety factor profile. 

After axisymmetric relaxation, three-dimensional modes are perturbed with random noise 
at very low amplitude. For this calculation we use a toroidal mesh with N# = 8 toroidal mesh 
points, so that three toroidal Fourier modes (n = 0, 1, 2) are included after dealiasing. This 
resolution is marginally acceptable for highly accurate calculations, but will demonstrate 
the utility of the TRIM algorithm for this problem. 

The kinetic energy in the n — 1 and n = 2 modes are shown in Fig. 41. The n = 2 mode 
is unstable and the n = 1 mode is stable, in agreement with linear calculations [29]. 

The linear eigenmode for the poloidal velocity is shown in Fig. 42, where the poloidal 
velocity vectors are shown at the eight toroidal locations included in the calculation. The 
flow pattern has the clear counterrotating vortex structure of an internal kink mode with 
dominant poloidal mode number m — 1 . This structure is seen to rotate twice around the 
torus, as required by an n = 2 mode. 

The instabilities computed here evolve on a fraction of the poloidal Alfven time, which is 
almost a factor of 10 longer than the toroidal Alfven time. Purely explicit methods require 
that the time step be taken at a fraction of the shortest time scale. In the example computed 
here we have used the semi-implicit method with a time step 30 times larger than that 
allowed by explicit numerical stability. 

6.4.3. Resistive Instability in a Torus 

The results presented in Sections 6.4. 1 and 6.4.2 described ideal instabilities, i.e., unstable 
normal modes of the ideal (infinitely conducting) MHD equations. However, some of the 
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V 

FIG. 39. The safety factor q(i/r) for the ITER equilibrium. 


most important instabilities exist only in the presence of finite resistivity. As a result of 
these modes the magnetic field configuration can change its topological properties, which 
are otherwise invariant in ideal MHD. These resistive instabilities [30] have no counterpart 
in ideal MHD and grow on a time scale that is a hybrid of the Alfven and resistive time 
scales. The study of these modes has been a primary focus of resistive MHD computations 
for over 20 years. 

We have begun to apply TRIM to resistive instabilities in a torus. We have studied the 
linear stability of a toroidal equilibrium whose stability properties are well known [31]. For 
this case, the poloidal cross section of the plasma is circular and the <7 -profile is given by 


q(r) = 2 


/ 1 + (r/r 0 ) 8 

v 1 + OW^o) 8 


1/4 


(84) 


where r is the minor radius (measured from the center of the circular poloidal cross 
section), r s2 is the radius of the q = 2 surface, and r 0 is the width of the current chan- 
nel. We use r s2 = 0.7 and r 0 = 0.6, which correspond to Run 1, Table 1 of Ref. [31]. The 
particular equilibrium has been supplied in numerical form [32]. In Ref. [31], the Lundquist 
number at the q = 2 surface was S = 2 x 10 4 . For our initial calculations, in order to 
save computer time, we have used an enhanced resistivity that gives S = 10 3 at the q = 2 
surface. 
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FIG. 40. Unstructured mesh for the ITER equilibrium. 


We perturb the initial equilibrium with random noise and consider a linear case with 
toroidal mode number n — 1 . We find an exponentially growing instability with a growth 
rate of ytAp =0.01 (where r a p is the poloidal Alfven time). This is to be compared with 
the result ytAp =0.017 of Ref. [31] with S = 2 x 10 4 . The effect of this resistive instability 
on the magnetic field topology is shown in Fig. 43, where we use the successive intersec- 
tions of four different magnetic field lines with the poloidal plane (a Poincare plot). The 
field line integration was performed with the code TUBE [33]. The magnetic islands that 
characterize resistive instabilities correspond topologically to a poloidal mode number 
m — 2, in agreement with Ref. [31]. 


7. SUMMARY AND DISCUSSION 

An algorithm for the solution of the time-dependent, primitive, resistive MHD equations 
in three-dimensional toroidal geometry has been developed. The algorithm uses an unstruc- 
tured, triangular mesh in the poloidal plane and a structured, pseudospectral method based 
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FIG. 41. Kinetic energy versus time for the n = 1 and n = 2 modes in the ITER equilibrium. 


on FFTs in the toroidal coordinate. This allows axisymmetric configurations with arbitrarily 
complex poloidal geometry to be accurately represented. Boundaries in the poloidal plane 
need not correspond to flux surfaces or be constrained by coordinate systems. The algorithm 
is fully conservative and maintains both the magnetic field and current density as solenoidal. 
Fluxes at cell interfaces are computed with a low-order upwind method. The semi-implicit 
method is used for time integration. 

The algorithm has been applied to four nonlinear test problems: a hydrodynamic shock 
tube; an MHD shock tube; supersonic flow past a sphere (both with and without an ini- 
tial dipole magnetic field); and, growth and saturation of toroidal instabilities. For both 
the hydrodynamic and MHD shock tube problems, good agreement with previous re- 
sults has been obtained. The primary inaccuracy is due to the numerical diffusion in- 
troduced by the low-order fluxes. Mesh adaption and refinement has been successfully 
applied to both the hydrodynamic and MHD cases and to the case of supersonic flow 
past a sphere. Linear growth and nonlinear saturation of three-dimensional kink modes 
in two analytic equilibria and in a highly elongated toroidal equilibrium representative 
of the ITER design have been computed. Time-step converged linear growth rates agree 
with previous linear stability calculations to within 20%. The nonlinear evolution of these 
modes has shown nonlinear mode coupling and spectral broadening and has demonstrated 
the utility of the semi-implicit method of time integration for these calculations. The lin- 
ear growth of a resistive tearing instability in toroidal geometry has also been calcula- 
ted. 

We briefly remark on the order of accuracy and efficiency of the algorithm. In Fig. 44 
we plot the numerical damping rate inherent in the algorithm as a function of number of 
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FIG. 44. Numerical damping rate vs number of triangles. 


triangles. These results were obtained during the relaxation of the Solov’ev equilibrium, 
as described in Section 6.4.1. We see that the numerical damping rate is approximately 
linear in 8 -- (A a s ) l/2 , confirming the first-order nature of the algorithm. This damping is 
consistent with the low growth rates found for these cases. The efficiency of the algorithm 
is displayed in Fig. 45, where we plot the CPU time required to compute for one e-folding 
time of a linear instability in the Solov’ev equilibrium as a function of At/Atcm., the ratio 
of the actual time step to the time step allowed by the explicit CFL condition. The behavior 
illustrated here is influenced by two factors. First, as the time step is increased above 
the CFL limit, the semi-implict coefficient, Eq. (72), becomes larger, and consequently 
more iteration are required to invert the semi-implict operator. Second, the e-folding time 
itself increases (the growth rate decreases) as the time step is increased (see Fig. 33). 
The most efficiency (i.e., the minimum CPU time per e-folding time) is obtained at about 
At/ A/cfl — 10. 

Several issues have arisen in the course of our investigation that we believe require further 
research. These are enumerated below. 

The first issue concerns the most efficient and accurate use of the primary and dual 
meshes, and the placement of dependent variables on them. In this work we have used the 
Delaunay triangles as the primary control volume and have chosen to define all volume 
densities (momentum, mass, and energy) at their centroids. All physical boundaries consist 
of triangle edges. We have not made use of the dual control volume elements consisting 
of the polygons centered at triangle vertices with edges connecting triangle centroids. (The 
edges of these polygons are used to define components of the magnetic field and cur- 
rent density.) With the present scheme differential operators such as the gradient and the 
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FIG. 45. CPU time per e-folding time for a linear instability in the Solov’ev equilibrium versus Ar/At CFL , 
the ratio of the time step to that allowed by the CFL condition. 


Laplacian acting on triangle-centered densities couple more than nearest neighbor triangles. 
Our experience with rectangular, structured meshes [16] indicates that the use of staggered, 
overlapping volume elements leads to the most compact, accurate, and physically motivated 
algorithms. In those methods, the pressure and momentum are not collocated as they are 
here, but are defined at the centers of the staggered primary and dual meshes, respectively. 
We have not experienced any severe problems that can be directly attributed to the non- 
compact formulation described in this paper, but the present algorithm seems to require 
somewhat more artificial viscosity to assure robustly stable computations than methods that 
use rectangular, structured meshes. Our initial attempts to formulate an exactly equivalent 
method using triangles and polygons has led to problems in consistently defining boundary 
edges and applying boundary conditions. 

The second issue concerns interpolation. In the present algorithm, interpolation from 
centroids to vertices, vertices to centroids, centroids to edges, edges to vertices, and vertices 
to edges are all required. The form of the interpolation can affect the accuracy and stability 
of the calculation. As described in Section 4.3, some of these interpolation schemes can 
become quite complicated and can lead to coupling beyond nearest neighbors. We have 
not devised a method for MHD using either rectangular or triangular meshes that does not 
require some interpolation (or averaging) from one grid to another. Interpolation is required 
during mesh refinement and coarsening. The number of interpolations per time step is also 
affected by the choice of primary and dual meshes. The issues of accuracy and required 
number of interpolations, and their affect on the performance of the algorithm, must be 
better understood. 
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The third issue concerns the criteria used for dynamic refinement and coarsening of the 
mesh. In Section 2.7 we presented criteria for adding or deleting triangles that are based on 
either the normalized average Laplacian operator or the normalized gradient. For the results 
presented in Section 6, we have used these criteria in conjunction with the mass density. 
Other criteria have been used in computational hydrodynamics [18]. Other problems may 
require still other criteria. For example, criteria based on current density may be useful 
for magnetic reconnection problems. The proper criteria for resistive MHD is yet to be 
determined. 

The fourth issue concerns the structure of the mesh itself. Delaunay triangles have many 
desirable properties. The algorithm is most accurate for equilateral triangles, for then the 
centroid dual mesh and the Voronoi dual mesh are equivalent. Deviations from equiangular- 
ity introduce errors, and highly obtuse triangles can lead to spikes in high order derivatives, 
such as current density. Mesh refinement can lead to large variations in triangle size and 
shape over the mesh, even though the Delaunay triangulation maximizes equiangularity 
in a global sense. These large variations can in turn affect the diagonal dominance of the 
Laplacian operator and cause the conjugate gradient algorithm to fail to converge. What 
is needed is a method for systematically redistributing the vertices in the poloidal plane 
to assure that all triangles are at least acute and that variations in mesh size are smoothly 
distributed in space. Clearly, more fundamental work needs to be done in this area. 

We have used low order approximations for the calculation of interface fluxes. The 
accurate computation of these fluxes has occupied the attention of computational 
hydrodynamicists for several years, and it has proven to be a crucial issue in the accu- 
rate engineering application of these methods. Improvements in these methods for MHD 
must eventually be addressed. 

Finally, we remark on the maximum Lundquist number (minimum dissipation) accessible 
by the TRIM algorithm. As is the case with all nonlinear calculations, we have found 
that, in the nonlinear regime, it is necessary to keep both the cell magnetic Reynolds’ 
number, Rma = v8/rj, and the cell viscous Reynolds’ number, R A = v8/v, of order unity. 
The Lundquist number is related to the cell magnetic Reynolds’ number by S = (v A /v) 
( L /8)Rm a • The maximum Lundquist number is therefore determined by both the particular 
problem under consideration (through the flow velocity and the Alfven velocity) and by the 
resolution (through the cell size 8, which is proportional to the square root of the number 
of triangles). Thus, for a case with v A /v ~ 10 and with 10 4 triangles (equivalent to a grid of 
100 x 100 on a rectangular mesh), S max ~ 10 3 . This is consistent with our experience. For a 
given problem, S max can only be increased by decreasing the size of the triangles. However, 
the scaling is not favorable: to increase S max by a factor of 10 with uniform triangle size, 
the number of triangles must be increased by a factor of 100 (^an order of magnitude for 
each spatial dimension). This is not to say that accurate but limited calculations cannot be 
performed at higher values of S. (See, for example, the ideal mode in the ITER equilibrium 
described in Section 6.4.2.) Indeed, the linear growth and initial nonlinear saturation of 
tearing modes in tokamaks can proceed for a time until the nonlinear coupling cascades 
significant energy to the shortest length scale <5. All that is required is that the linear mode 
structure near the resonant surface be well resolved, which can be accomplished by packing 
grid cells there. However, if Rma > T this energy will eventually and inevitably build to 
a level sufficient to render that calculation meaningless. This is the case for all nonlinear 
algorithms, not just the one described here. 
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particular, evolution of equilibrium NC, which maintains 
its large-scale magnetic structure and evolves on the fastest 
timescale, may provide a viable physical model for compact 
loop flares. FC models, which produce weaker current 
sheets and dissipate smaller amounts of energy, and evolve 
on longer timescales, would be expected to contribute to 
coronal heating, and may explain the long-lived hot loops 
observed by Yohkoh. 

In summary, line tying has two contrasting effects on the 
behavior of coronal loops. For sufficiently short loops, it 
acts as a stabilizing agent to prevent perturbations from 
growing, allowing loops to build up and store magnetic 
energy. Once the loops become unstable, line tying can 
cause current sheets to form, resulting in an enhanced level 
of energy dissipation compared to non-line-tied configu- 
rations. It would be worthwhile to extend the calculations 
presented here to include a more realistic description of the 
thermodynamics, and to include the effect of toroidal curva- 
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ture, which may cause a rapid outward expansion of loops 
as they are twisted (Amari et al. 1996). 
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of the very rapid increase of growth rate as marginality is 
crossed. 

4. NONLINEAR RESULTS 

We have performed three simulations to study the non- 
linear evolution of kink instabilities that develop in the FC 
and NC fields. The MHD equations (l)-(4) are advanced 
numerically in time by means of a three-dimensional semi- 
implicit algorithm, which makes use of a Fourier expansion 
in the poloidal direction and of a finite-difference technique 
in the axial and radial directions. We define a radial mesh 
extending between r = 0 and r = R, where we set a conduct- 
ing boundary. Since we expect to find current concentra- 
tions in the vicinity of the center of the cylinder, the mesh is 
designed as in Mikic et al. (1990), with greater resolution for 
small r’s and increasingly sparser at greater radii: 13% of 
the points are within r = 0.05 R, 26% are within r = 0.1R, 
55% are within r = 0.25 R, and 82% are within r = 0.5 R. 
The ratio of the largest Ar to the smallest is 13.9. R must be 
large enough to avoid artificial stabilizing effects on the 
evolution of the system. We choose R = 20 for the FC 
fields, which develop global kink modes, and R = 12.5 for 
the NC equilibrium, in which the internal kink mode is 
effectively confined by the outer constant axial field. The 
axial and poloidal meshes are taken to be uniform in both 
cases. A 63 x 63 x 32r — z — 6 mesh is employed. We use a 
dealiased Fast Fourier Transform algorithm, so that only 
modes with 0 < m < 10 are retained (equivalent to § of 
available Fourier space). The inversion of self-adjoint 
matrices, generated by the semi-implicit formulation and by 
the fully implicit treatment of the viscous and resistive 
terms, is accomplished through a conjugate gradient solver 
with diagonal preconditioning. The details of the code are 
described elsewhere (Lionello, Mikic, & Schnack 1998). The 
initial conditions assume constant mass density, p = 1, con- 
stant temperature, a plasma beta = 0.01, which is typical 
of the solar corona and an m = 1 perturbation of amplitude 
v ~ 6 x 10 _ 4 is applied to the three equilibria described 
above. In order to save computing time the radial and axial 
profiles of the initial perturbation are derived from the 
linear code. 

The computation is started with no resistivity and with a 
small value of the viscosity, v ~ 10 -3 , to be able to follow 
the entire ideal linear phase of the instability. We find that 
in all three cases during the first stage of the nonlinear ideal 
evolution small scales are created in all three directions and 
eventually resistivity must be introduced in the computation 
in order to dissipate the energy which accumulates at small 
scales and to avoid the development of numerical insta- 
bilities (Einaudi et al. 1997). Therefore after a time T r , whose 
value depends on the equilibrium under investigation, the 
resistive terms are turned on with rj = 10~ 2 ( T r = 195, 160, 
and 115, respectively, for FC1, FC2, and NC), and at the 
same time the viscosity is set to the same value. Such a value 
is necessary to describe the behavior of small scales in our 
moderate resolution simulations properly. Resistivity is 
actually important only locally where small scales are 
created by the nonlinear dynamics. Attempts to link the 
inclusion of resistivity to the magnitude of the currents (to 
avoid the diffusion of the large-scale structures) in our code 
were unsuccessful, probably because of the implicit treat- 
ment of the dissipative terms. 

Once resistivity is turned on, the instability eventually 
saturates in a timescale that appears to be independent of 


the equilibrium adopted even though the linear timescales 
are different. Such a behavior is well described by the time 
history of the magnetic and kinetic energies contained in the 
different poloidal modes, which is shown in Figure 6 (note 
the cusps in the profiles at the time when resistivity is turned 
on). 

In all runs an m = 0 transient mode shows up in the 
kinetic energy history, which is due to adaptation of the 
analytical equilibria to the grid. We have derived the 
approximate growth rate of the m= 1 mode in the three 
cases, finding a satisfactory agreement with the growth rates 
obtained using the linear code presented in Figure 3. From 
Figure 6 it is also clear that the growth of the higher m 
modes is due to forcing from the m = 1 mode via quadratic 
nonlinear terms, since, e.g., the m = 2 mode, once the m = 1 
mode is well into the exponential growth phase, also 
exhibits exponential growth with a rate that is twice as large 
as that of the m = 1. 

In order to investigate the distribution of the currents in 
the loop and to better understand the physical meaning of 
the creation of small scales, we have computed the modulus 
of the current density at different heights for the three equi- 
libria at the beginning of the resistive phase, and we present 
the results for z = 0 and z = f L in Figure 7. It is evident 
that the introduction of the resistivity is necessary to cor- 
rectly describe the evolution of current peaks that have been 
formed during the ideal nonlinear evolution. The overall 
topology of the current structure differs considerably from 
one equilibrium to the other, reflecting the different physical 
nature of the current concentration. We have already 
described how for model FC1 the kink extends far into the 
regions of the original potential field and the current dis- 
tribution does not show any strong concentration in the 
linear phase. In the nonlinear phase two current peaks 
develop, one per each half loop, the location of the peaks 
depending on the total length of the loop. This fact can be 
seen comparing the results for FC1 with the results present- 
ed in Velli et al. (1996) for their model a, which has the same 
radial structure as FC1 and length L = 25. The current 
concentrations at z = 0 are just the “ tips ” of the two layers 
developed in the legs of the loop. The formation of such 
current layers is a direct consequence of the fact that the 
magnetic lines cannot move at photospheric levels produc- 
ing strong gradients where the kink folds on itself in each 
leg of the loop. 

The presence of “ resonant ” regions in the linear phase 
can be another possible cause for the formation of current 
concentrations in the subsequent nonlinear phase. In a 
“resonant” region the ideal terms determining the time 
evolution of the radial component of the perturbed mag- 
netic field go through zero, generating a current concentra- 
tion within the region. For both FC2 and NC models the 
“ resonant ” region is located at the loop apex and current 
sheets develop there. The difference between the two models 
is that in the NC case only the internal part of the plasma 
column moves outward, piling up magnetic energy at the 
border with the potential region and contributing to the 
local enhancement of the current. In the FC2 case the radial 
displacement is more extended and therefore the resonant 
region is carried by the flows outward with a smaller oppo- 
sition of the potential field. 

In all cases current layers are formed and the evolution 
becomes resistive with a consequent occurrence of magnetic 
reconnection (see also Baty & Heyvaerts 1996). 
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Fig. 6. — Poloidal mode histories of the magnetic and kinetic energy for equilibria FC1 ( top panel), FC2 ( middle panel), and NC ( bottom panel) 


In Figures 8 and 9 we present field line plots for the FC 
and NC models at the beginning of the resistive phase and 
at saturation of the instability when all poloidal modes start 
declining. The lines start from the bottom of the loop at 
different radii (r = 2, 4) and angles, in the region where the 
current flows, and are integrated until they reach the top. 
This choice of field line imaging makes the representation 
asymmetrical, though clearly a perfectly symmetrical image 
would appear if we plotted also the corresponding field lines 
originating from the top of the loop. The axial dimensions 
of the structures are proportional to the real lengths of the 
loops. Each magnetic line connects two fluid elements that 
are rooted in the photosphere at opposite sides of the loop. 
During the ideal phase the field line can only be bent and 
twisted since the fluid elements in the photosphere cannot 
move. In the resistive phase the field lines can reconnect, 
altering the field topology and therefore the connectivity 
between fluid elements in the photosphere. 

The fact that reconnection has occurred during the 
resistive nonlinear evolution of all three models is evident 
from a comparison of left and right panels in Figures 8 and 
9. Following the previous discussion on the properties of 
the dynamics of the three models, it is not surprising to 
notice that the final magnetic topologies of the three models 
differ considerably even if the physical mechanism deter- 
mining the evolution is the same in all cases. This is due to 
the fact that the current layers resulting from the ideal non- 
linear evolution are located in different regions of the loop. 


Reconnection starts influencing those magnetic lines that go 
through the current layer and therefore changes the connec- 
tivity of those fluid elements at the photosphere where such 
lines originate. Away from the current layer the lines are 
frozen and move with the local velocity of the fluid. 

In the case of model FC1 magnetic lines are affected by 
reconnection in both legs of the loop and the kink progres- 
sively involves field lines at greater and greater distances 
from the axis of the loop. As a result lines with an originally 
small amount of twist are connected with lines originally in 
the potential region of the loop with a much greater twist 
and eventually the original magnetic field topology is com- 
pletely destroyed, leading to what appears to be an ergodic 
field. We have not checked whether there are field lines that 
are actually space filling, or field lines that are totally dis- 
connected from the photosphere, but an indication may be 
given by the behavior of the function a =j • B/B 2 , which is 
a constant along field lines. Initially, this quantity is non- 
vanishing only within the current channel and monotoni- 
cally decreasing outward (it does not depend on z). At 
saturation, this function appears to develop a plateau over 
the full extent of the perturbed region. At the same time, the 
field lines are clearly not the laminar constant a force-free 
field. 

Models FC2 and NC develop a current layer at the loop 
apex and suffer a less extended kink instability. In the NC 
case the kink is totally confined within the current channel 
and the original potential field contains no twist. For model 
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Solanki et al. (1993). According to Hagyard et al. (1983) and 
others, the observed vertical field gradient ranges from 0.2 
G km “ 1 to several G km “ 1 (for reference see Balthasar & 
Schmidt 1993). Hagyard et al. (1983) suggest 0.2 G km -1 , 
while Balthasar & Schmidt (1993) suggest 2-3 G km 1 and 
Pahlke & Wiehr (1990), 2 G km -1 . Balthasar & Schmidt 
(1993) obtained the nearly independent vertical field gra- 
dient on location within a sunspot. In the present study, we 
have considered four different cases: dB z /dZ = 0 G km -1 , 
0.5 G km -1 , 1 G km -1 , and 2 G km -1 . By virtue of equa- 
tions (9) and (11), the longitudinal field configuration refer- 
enced to a common geometrical depth (Z = Z D ) can be 
obtained by 

£ Z (R, Z D ) = B z (a) + AZ(oc)(dB z /dZ) . (12) 

Here we express our new shape function as 

D(cc) = D( 0) exp ( — a n ) = D( 0) exp [-(a p R/R p Tl , (13) 

where n and ol p are free parameters that will be chosen to fit 
D(oc)/D(0) to the derived B Z (R, Z D )/B z (0 , Z D ) through 
equation (2). The old shape function is a Gaussian (n = 2) 
with ol p = 1.63, which was inferred from the longitudinal 
field configuration referenced to optical depth. Since we 
have converted the magnetic field distribution referenced to 
optical depth to one referenced to geometrical depth, we 
have to derive new shape functions that will conform to the 
new prescription. For this work we generalize the shape 
function to that given by equation (13). This corresponds to 
an improvement of the upper boundary condition of Yun’s 
sunspot models. 

A set of the two numerical parameters n and cc p has been 
determined by means of nonlinear least-squares fitting for 
each of the four different cases considered above. The 
resulting values are tabulated in the third and fourth 
columns of Table 1. In the table we may note that case 1 
(zero field gradient) is very similar to the old one chosen by 
Yun (1968). The total flux O and the parameter / conform- 
ing to the newly derived shape function have been estimated 
by numerically solving equations (3) and (5). We define A 
and C as A = f/D(0) and C = ®/nD(0) and list the resulting 
values of A and C for the four cases in Table 1. Finally, the 
radial component of the field B r (R , Z D ) is computed from 

B r (R 9 Z d ) = aD(am/dZ) ZD 

= exp [ — (ai?/R p )"](tan ^ p )(a/a p ) . 


2.3. Boundary Condition 

The upper boundary condition is a free parameter, which 
can be specified by the inclination angle \j/ p at the outer edge 
of the penumbra at a given geometrical depth Z D . The 


upper boundary condition is given by 


y\z D ) 


y\Z D ) tan ij/ p 

[D(0)] 1/2 a p 


(15) 


where Z D is taken as the depth of t = § at the sunspot 
center and x l/ p is the angle of inclination of a line of force 
passing through the outer edge of the penumbra at the 
geometric depth Z D . The value of i \/ p here should be smaller 
than the observed angle of inclination at the outer edge of 
the penumbra. However, Yun (1970) neglected the Wilson 
effect and made use of values based on the observations. 
The lApmax * n Table 1 is the upper limit of \ j/ p9 below which 
no off-centered local maximum of the total field strength is 
found. For all cases, the values of ^ pmax are smaller than the 
observed inclination angle ^ obs (0). 


3. RESULT AND DISCUSSION 

With the new shape functions derived in the previous 
section (see Table 1), four sets of magnetostatic sunspot 
models with T spot = 4000 K and 4400 K have been con- 
structed, in each of which i jf p is varied to fit the observ- 
ations. In Figure 2, the field strengths at the spot axis 
computed for three values of i jf p are compared with the 
observed ones ( solid line ) determined by the empirical rela- 
tion (eq. [6]). As can be seen from the figure, the computed 
field strengths match well the observed one for a particular 
x j/ p , which is named as the most desirable one. We tabulate 
the most desirable i \f p in Table 1 for comparison with the 
upper limit angle ^ pmax . As can be seen in Table 1, the most 
desirable ones of all cases, except for case 1 (zero field 
gradient), do not exceed the upper limit i A pmax , beyond 
which the off-centered maximum field strength begins to 
show up. It implies that the shortcomings of Yun’s model 
do not appear in the computed models since the upper 
boundary condition is different from that of Yun (1970) by 
taking into account the Wilson effect. In addition, the deter- 
mination of the most desirable value of i j/ p for a assumed 
field gradient results in the models that are essentially char- 
acterized only by the effective temperature T spot at the spot 
axis. This makes it possible to investigate the dependence of 
the physical characteristics of model spots on T spot and/or 
on the spot size. The physical characteristics of the com- 
puted model sunspots are summarized in Table 2, which are 
specified by the most desirable t ]/ p — 50°, 63°, 67°, 71°, and 
T spot = 4000 K, 4400 K, respectively. As can be seen in 
Table 2, thermodynamic parameters P g and p depend only 
on the thermal structures of model spots. 

So far, the vertical field gradient has been assumed to be 
independent of the position within a sunspot. In order to 
examine the radial dependence effect of the magnetic field 


TABLE 1 

Model Parameters Estimated for Different Vertical Field Gradients 


P <Aobs 

dB ? JdZ lApmax (Z D ) (0) 

Model (G km -1 ) n a p A C (deg) (deg) (deg) 


Yun 0 2.00 1.63 0.50 1.00 66.5 

Case 1 0 1.86 1.63 0.49 1.03 69.1 71 72 

Case 2 0.50 1.65 1.43 0.48 1.11 68.7 67 72 

Case 3 1.00 1.49 1.22 0.47 1.19 66.8 63 72 

Case 4 2.00 1.24 0.81 0.47 1.44 58.6 50 72 

Case 5 — 0.5a + 1 1.70 1.39 0.48 1.09 67.6 67 72 
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Fig. 2. — Computed maximum field strengths (points ) at the photosphere of sunspots vs. the observed ones ( solid line) given by eq. (6). dB JdZ = (a) 0 G 
km -1 , (b) 0.5 G km -1 , (c) 1 Gkm -1 , and (d) 2 Gkm -1 . The two field strengths are consistent with each other for a particular ^ Three i k for each case are 
adopted for comparison. 


gradient, we have taken a case characterized by 

dB z /dZ = -0.5a + 1 G km" 1 , (16) 

where we assumed that the vertical field gradient decreases 
with the radial distance from the spot center. With a new 
shape function obtained from the field configuration char- 
acterized by equation (16), we have calculated the sunspot 
models by the same method as in the above four cases and 
included the results in the last rows of Tables 1 and 2. As 
seen from the tables, the physical characteristics of the com- 
puted models in case 5 are nearly similar to those of case 2, 


where the vertical field gradient was taken as 0.5 G km -1 . 
This implies that ij/ p depends critically on the vertical gra- 
dient near the outer part of the penumbra. 

Figure 3 shows the field distribution of the radial, verti- 
cal, and total magnetic field strengths at the depth of Z d 
calculated by using new shape functions and equation (14) 
for four different field gradients (0.5, 1, 2, —0.5a + 1 G 
km -1 ). Since the inclination angle of case 4 ( Fig. 3c) 
changes too much from 50° at Z D to 72° at Z = 0, it could 
not be compatible with the observation of deep sunspot 
penumbra (Solanki & Schmidt 1993). Among the con- 


TABLE 2 

Physical Characteristics at the Center of Model Sunspots at the Surface Z d 


T <D B 

Model (K) (deg) (Mx) (G) 


Case 1 4000 71 2.4E22 3337 

4400 71 5.6E21 2720 

Case 2 4000 67 2.4E22 3406 

4400 67 5.6E21 2774 

Case 3 4000 63 2.4E22 3387 

4400 63 5.6E21 2759 

Case 4 4000 50 2.4E22 3377 

4400 50 5.6E21 2759 

Case 5 4000 67 2.4E22 3374 


4400 67 5.6E21 2748 


p g 

(dynes cm -2 ) 

p 

(gem 3 ) 

(dB/dZ) c 
(G km _ *) 

(Z D ) 

(arcsec) 

1.99E5 

7.86E-7 

0.79 

33.5 

1.73E5 

6.22E-7 

1.20 

17.9 

1.99E5 

7.86E-7 

0.79 

28.0 

1.73E5 

6.22E-7 

1.20 

15.0 

1.99E5 

7.86E-7 

0.79 

23.2 

1.73E5 

6.22E-7 

1.20 

12.4 

1.99E5 

7.86E-7 

0.79 

14.0 

1.73E5 

6.22E-7 

1.20 

7.5 

1.99E5 

7.86E-7 

0.79 

27.6 

1.73E5 

6.22E-7 

1.20 

14.8 
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Fig. 3. — Magnetic field configuration characterized by new shape functions and the most desirable t j/ p (see Table 1). The most desirable ij/ p was determined 
by comparing the computed field strengths with the observed ones. dB z /dZ = (a) 0.5 G km - 1 , ( b ) 1 G km - 1 , (c) 2 G km~ \ and ( d) —0.5a + 1 G km - 1 . 


sidered models, cases 3 and 5 seem to be the most self- 
consistent model in the sense that the resulting (dB/dZ) c 
(0.8-1.2 G km -1 ) comes out closer to the assumed vertical 
field gradient of 1 G km -1 . Considering the recent obser- 
vation by Bruls et al. (1995), who found that the vertical 
gradient of magnetic fields in the penumbra slowly declines 
as a function of radial distance, we suggest that case 5 
should be a most reasonable one. It is noted that the com- 
puted d\l/ p /dZ ranges from 0?01 km -1 to 0?03 km -1 at the 
outer edge of the penumbra, which falls within the range 
suggested by Solanki et al. (1993). 

Finally, we note that in our model the amount of the 
Wilson depression appears to have a linear relationship 
with T spot /T ph and T spot /T ph can be related to R u , a directly 
observable quantity with the use of equations (6) and (7). In 
Figure 4 we plot the Wilson depression calculated from 
models as a function of R u to get their linear relationship as 

WD(km) = 26.4£ M [arcsec] + 356.2 (17) 

for R u > 7". We suggest this model-predicted relation could 
be tested by high-resolution observation in the future. The 
amount of the Wilson depression we presented in Figure 4 
lies within the range of those suggested by Solanki et al. 
(1993) and Gokhale & Zwaan (1972). 

The self-similarity assumption by ST58 is adopted in our 
study to make the problem mathematically tractable. There 



Fig. 4. — Computed Wilson depression as a function of the umbral size. 
Crosses are the computed values, and the solid line is their linear regres- 
sion. 
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is no known physical reason for the sunspot field configu- 
ration to be self-similar. However, we can at least assert that 
the self-similarity, if any, can hold only in the region of high 
plasma /?, where the plasma pressure can confine the mag- 
netic field. So far no direct measurement of the sub- 
photospheric field configuration has been made. Lee et al. 
(1993) investigated the radial magnetic field structure of a 
single isolated sunspot at the coronal base using microwave 
spectroscopy. They showed that the radial field distribution, 
at least inside the inner penumbra, has a Gaussian form and 
is very similar to that in the photosphere. If the field con- 
figuration is more or less self-similar even above the photo- 
sphere, we may expect that it would be more so in the 
subphotospheric region. 

In the present study, we have not made a closer exami- 
nation of the optical properties of the models. A continuum 
radiative transfer analysis of the emergent flux at various 
center-to-limb angles would give us a more definitive under- 
standing of the Wilson effect. However, it is beyond the 
scope of this paper and is reserved for future studies. 

The sunspot models in this study are constructed under 
the assumption of the azimuthal symmetry. Title et al. 
(1993) found from their high-resolution observations at the 
Swedish Solar Observatory in La Palma that the inclina- 
tion angle of the magnetic field in the penumbra of sunspots 
oscillates rapidly with azimuth. In relation to this, Martens 
et al. (1996) presented a constant-a force-free model for the 
magnetic field in fluted sunspots. Our symmetric models do 
not accommodate any azimuthal variations, but are mainly 
concerned with the mean-field structure of isolated sun- 
spots, especially in the high-/? region in and below the 
photosphere. 

4. SUMMARY AND CONCLUSION 

In the present study, we have revisited Yun’s sunspot 
models by taking into account the effect of the Wilson 
depression. For this purpose, we represented the observed 
field distribution of sunspots with the Skumanich dipole 
model and converted it to one referenced to a common 
geometric depth. In converting the Skumanich dipole 
model to the geometrically referenced one, we have made 
use of the radial dependence of the Wilson depression sug- 
gested by Solanki et al. (1993). With the use of the Wilson 
depression of Solanki et al. (1993) along with five different 
vertical gradients (0, 0.5, 1.0, 2.0, and — 0.5a + 1 G km - x ), a 
set of new shape functions has been derived by fitting them 
to the geometrically referenced longitudinal field distribu- 
tions. The shape functions are assumed to be represented by 
D( a) = D( 0) exp ( — a”). To reduce the number of free param- 
eters in the model computations, we have utilized a recent 
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empirical relation between maximum field strength and 
effective temperature of spots given by Kopp & Rabin 
(1992), together with the relation between the field strength 
and the umbral size suggested by Kopp & Rabin (1992) and 
Martinez Pillet & Vazquez (1993). With the new shape func- 
tions conforming to geometrically referenced field configu- 
rations characterized by the five different vertical field 
gradients, five sets of magnetostatic sunspot models have 
been computed by varying t j/ p . 

The main result in this study is summarized as follows : 

1. The computed models for a given field gradient can be 
characterized only by 7^ pot , the effective temperature at the 
center of spots. The most desirable values of \j/ p , the upper 
boundary condition, for the five different cases were deter- 
mined by comparing the computed maximum surface field 
strengths with the observed ones. 

2. Especially, the computed models do not have the 
shortcoming raised by Osherovich (1982), demonstrating 
that the shortcoming is not inherent to the self-similar 
models. The similarity assumption would break down in the 
higher solar atmosphere due to the surrounding weak gas 
pressure as commented by Steiner, Pneuman, & Stenflo 
(1986) and Pneuman, Solanki, & Stenflo (1986). 

3. The resulting sunspot models support the observed 
empirical relations given by Kopp & Rabin (1992) and 
Martinez Pillet & Vazquez (1993). 

4. The computed d\j/ p /dZ ranges from 0?01 to 0?03 km - 1 
at the outer edge of the penumbra, which falls within the 
range suggested by Solanki et al. (1993). 

5. The present study supports the self-similarity for the 
lateral magnetic vector structure of fairly isolated sunspots 
found by Keppens & Martinez Pillet (1996). 

We presented in Table 2 the physical characteristics of T, 
i j/ p , <I>, B , P g , p, (dB/dZ) c , and R p (Z D ) expected in each 
model. It is concluded that when the Wilson depression is 
properly taken into account, the similarity field configu- 
ration represents quite well the structure of sub- 
photospheric sunspots. 
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